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Preface 


"The  Second  U.S. -Japan  Workshop  on  Advanced  Plasma  Modeling" 
meet  at  IPP  Nagoya  from  March  23-27th,  1987.  The  program  for  the 
meeting  is  attached  as  appendix  I. 

There  were  seven  American  participants  and  26  Japanese 
participants ; the  list  of  participants  is  attached  as  appendix  II. 

We  consider  that  this  workshop  was  very  successful.  There 
was  lively  discussion  from  both  sides  on  detail  of  modeling 
techniques,  their  advantages  and  disadvantages,  how  they  might  be 
improved  and  what  types  of  physics  problems  could  be  treated.  It 
was  clear  that  Japanese-U . S . collaboration  has  played  a  large  rcle 
in  developing  many  powerful  new  techniques  for  handling  large 
time  step  large  space  scale  problems.  There  were  also 
developments  that  took  place  in  one  country  or  the  other  and  it 
was  clear  that  the  participants  learned  a  lot  from  hearing  of 
each  others  work.  We  have  drawn  a  number  of  conclusion  from  this 
workshop . 


J.M.  Dawson 


T.  Kamimura 


Summary 


Gyrokinetic  Particle  Simulation 

The  numerical  properties  of  a  gyrokinetic  plasma  have  been  discussed  in  detail 
in  this  workshop.  It  is  generally  agreed  that  considerable  improvement  in  time 
step,  grid  spacing  and  noise  level  can  indeed  be  realized  with  the  present  model. 
Numerical  schemes  for  solving  the  nonlinear  gyrokinetic  equations  have  also  been 
presented.  The  consensus  is  that  the  treatment  of  the  electron  motion  parallel  to 
the  magnetic  field  is  very  crucial.  The  proposed  sub-cycling  scheme  should  be 
examined  carefully. 

With  the  arrival  of  the  present  generation  of  super-cornputerSyProf.  Dawson 

- '  V  y  f  c  e> 

feels  very  strongly  that  it  is=the  time  fer-n^to  work  actively  toward  the  realistic 
simulation  of  tokamak  discharges  with  particle  codes  —  gyrokinetic  or  implicit. 
However,  there  are  numerous  important  issues  one  has  to  resolve  before  achieving 
that  goal,  such  as  particle  re-cycling  techniques,  3D  diagnostics,  the  MHD 
equilibrium  problem,  the  development  of  collision  operators  and  the  choice  of  the 
appropriate  coordinate  systems  for  the  simulation.  Another  important  aspect  is 
the  inclusion  of  the  ^kitchen  ph^sics^  in  the  code.  _ 

In  view  of  the  man  power  shortage  in  this  area,  Profs.  Dawson  and  Kamimura 
have  proposed  a  collaboration  effort  between  the  US  and  Japan.  The  initial  phase 
of  this  joint  endeavor  will  start  this  summer  at  UCLA. 

W.  W.  Lee 


-  2  - 

'  '  •  -*:s* 


A 


Append? 


x  I 


US-Japar  workshop  on  Advanced  Plasma  Mode  1 i ng 
March  23-27,  1987 
IPP,  Nagoya,  Japan 


AGENDA 

March  23  (Monday) 

Registration  (9:30) 

Opening  Remarks  (10:00-10:20) 

J.  M.  Dawson  and  T. Kamimura 

Morning  Session  (10:20-12:30)  Chairman  I.Kawakami 

A. Aydemir(IFS)  :  Generalized  Reduced  MHD  and  Full  MHD  Calculations 

with  Semi-Implicit  Techniques 

D. D. Schnack(SAI)  :  Semi-Implicit  Methods  for  3D  MHD  Computations 
G. Kurita(JAERI)  :  Nonlinear  Evolution  of  Free  Boundary  Modes  in  Tokamak 

Lunch  (12:30-14:00) 

Afternoon  Session  ©  (14:00-15:30)  Chairman;  A. Aydemir 

T. Hayas  hi  (H I  FT)  :  3D-Equi  1  ibrium  and  Stability  code  of  Helical  System 

T. Hayas hi (HIFT)  :  A  Simulation  of  Driven  Reconnection 

by  a  high  precision  MHD  code 

Coffee  Break  (15:30-15:45) 

Afternoon  Session  ®  (15:45-17:15)  Chairman:  G.  Kurita 

T.Ogino  :  An  MHD  Simulation  OP  the  Solar  Wind  and  Gome t a ry  Plasma 

(Nagoya  Univ.)  (An  MHD  Model  with  Plasma  Production) 

C. Z. Cheng(PPPL)  :  NOVA-2:  A  Kinetic  MHD  Stability  Code 
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March  24  (Tuesday) 


Horning  Session  ©  (9:15-10:50)  Chairman:  H.  Abe 

J.  H. Dawson(UCLA)  :  A  Hybrid  Vlasov-Fluid  Hodel  with  Kinetic  Ions 

and  Massless  Fluid  Electrons 

W. W.  Lee (PPPL)  ;  Gyrokinetic  Particle  Simulation  of  Finite-Beta  Plasma 

Coffee  Break  (10:50-11:00) 

Horning  Session  0  (11:00-12:30)  Chairman:  C.Z. Cheng 

H. Tanaka(HIFT)  :  Macroscale  Particle  Simulation  of  Kinetic  Alfven  Waves 

(  Code  Description  and  Application  ) 

A. Friedman(LLL)  :  Advanced  Particle-in-cell  Plasma  Modeling  at  LLNL 

Lunch  (12:30-14:00) 

Afternoon  Session  ©  (14:00-15:50)  Chairman:  M. Tanaka 

K. Nishihara  :  Particle-Particle  Particle-Mesh  code 

(Osaka  Univ. )  for  Nonideal  High  Density  Plasma  and  Its  Application 

on  Rayl 8 i gh-Tay lor  Instability  in  ICF  Plasma 

S.Kawata  :  TRIPIC :  Triangular-Mesh  PIC  code 

(Nagaoka  Tec.)  for  LIB  Diode  Simulation 

S.Y.Kim  :  Particle  Simulation  Technique  using  FEM  Methods 

(KIT,  Korea) 

Coffee  Break  (15:50-16:00) 

Afternoon  Session  ®  (16:00-17:30)  Chairman:  M.Aizawa 

K.Miki  :  A  Domain  Decomposition  and  Overlapping  Method 

(Hitachi  Ltd.)  for  3D  Large  Scale  Numerical  Simulation 

R. Sydora(UCLA)  :  Subtraction  Technique  for  Plasma  Physics  and 

Evaluation  of  Numerical  Effects  in  Codes 

Group  Dinner  (18:30-  ) 
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March  25  (Wednesday) 


Horning  Session  (9:15-10:50)  Chairman:  T.Ogino 

A.  Friedman(LLL)  :  Numerically  Induced  Stochasticity 

T. Yabe (Osaka  Univ.)  :  Reduced  Atomic  Model  for  Calculation  of  Charge 

Distribution  in  Multiply-Charged  Plasma 

Working  lime  for  Discussion  Leaders  and  Group  Photo  (10:50-11:30) 

Lunch  (11:30-13:00) 


Discussion  and  Working  Sessions  ©  (13:00-17:30) 

©  Advanced  MHD  Model  (13:00-15:00) 

Discussion  Leaders:  T.  Hayashi  and  D.  D.  Schnack 

comments:  I.Kawakami,  T.Ogino,  and  T. Hayashi 

Coffee  Break  (15:00-15:30) 

©  Implicit  Particle  Model  (15:30-17:30) 

Discussion  Leaders:  K. N is  hi hara  and  A. Friedman 

comments:  T.Kamimura,  H.Abe,  M.  Tanaka,  A.  Friedman,  and  W.  W.  Lee 

March  26  (Thursday) 

Discussion  and  Working  Sessions  ©  (9:15-15:00) 

©  Gyrokir.etic  Particle  Model  (9:15-11:15) 

Discussion  Leaders:  H.Naitou  and  W.  W.  Lee 

comments;  W. W. Lee  and  J. M. Dawson 

Coffee  Break  (11:15-11:30) 
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@-l  Improvement  of  Simulation  Model  (  1  1:30  -  1  2:30  ) 

Discussion  Leaders:  T. Yabe,  I. Katanunia,  and  R.Sydora 

comments;  S.Takeuchi,  K.Hanatani,  and  V. Abe 

Lunch  (12:30-14:00) 


@-2  Continuation  (  14:00-1  5:00) 


Summary  and  Conclusions  (15:00-16:00) 


Co-chairmen  :  J . M. Dawson  and  T.  Kamimura 


March  21  (Friday) 


Working  Session 


Closing 


The  maximum  time  allowed  for  speakers  including  discussions; 
(40  -  45)  min.  for  the  oral  presentations 
(10  -  15)  min.  for  comments 


Appendix  II 


List  of  Participants 

US-Japan  Workshop  on  Advanced  Plasma  modeling 
March  23-27,  1987 
IPP,  Nagoya,  Japan 


U.S.A. 

Dawson  John  M.  University  of  California,  Los  Angeles 

Department  of  Physics  405  Hilgard  Ave. 
Los  Angeles,  Calif.  90024 

Aydemir  Ahmet  University  of  Texas  Austin 

Institute  for  Fusion  Studies 
Austin,  Texas  78712 

Cheng  C.Z.  Plasma  Physics  Laboratory 

Princeton  University 

P.O.  Box  451,  Princeton,  N.J.  08544 

Friedman  Alex  University  of  California 

Lawrence  Livermore  Lab. 

P.O.  Box  808,  Livermore  CA.  94550 

Lee  W.W.  Plasma  Physics  Laboratory 

Princeton  University 

P.O.  Box  451,  Princeton,  N.J.  08544 

Schnack  Dal ton. D. 

Science  Applications  International  Corp. 
San  Diego,  Calif.  92121 

Sydora  Richard  University  of  California,  Los  Angeles 

Department  of  Physics 
405  Hilgard  Ave. 

Los  Angeles,  Calif.  90024 
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JAPAN 


A 


Abe  Hirotada  Dept,  of  Engineering 

Kyoto  Univ.,  Kyoto  606 

Abe  Yoshihiko  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Aizawa  Masamitsu  Atomic  Energy  Research  Institute 

College  of  Science  &  Technology 
Nihon  Univ.,  Tokyo  101 

Amano  Tsuneo  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Hanatani  Kiyoshi  Plasma  Physics  Laboratory 

Kyoto  Univ.,  Gokasho,  Uji,  Kyoto  611 

Hayashi  Takaya  Institute  for  Fusion  Theory 

Hiroshima  Univ.,  Hiroshima  730 

Ichiguchi  Katsuji  Plasma  Physics  Laboratory 

Kyoto  Univ.,  Gokasho,  Uji,  Kyoto  611 

Ichikawa  Yoshi  H.  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Ishiguro  Seiji  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Kamimura  Tetsuo  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Katanuma  Isao  Institute  of  Applied  Physics 

Plasma  Research  Center 
Univ.  of  Tsukuba,  Ibaraki  305 

Kawakami  Ichiro  College  of  Science  £  Technology 

Nihon  Univ.,  Tokyo  101 
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Kawata  Shigeo  Technological  Univ.  of  Nagaoka 

Nagaoka,  Niigata,  940-21 

Kurita  Genichi  Japan  Atomic  Energy  Research  Institue 

Tokai  Research  Establishment 
Ibaraki  319-11 

Miki  Kazuyoshi  Energy  Research  Laboratory 

Hitachi,  Ltd., 

Hitachi-shi,  Ibaraki,  316 

Mizuno  Yukio  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Naitou  Hiroshi  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  464 

Nakajima  Noriyoshi 

Institute  of  Plasma  Physics 
Nagoya  Univ.,  Nagoya  464 

Nishihara  Katsunobu 

Institute  of  Laser  Engineering 

Osaka  Univ.,  Suita  565 

Ogino  Tatsuki  Research  Institute  of  Atmospherics 

Nagoya  Univ.,  Toyokawa  442 

Takemoto  Yukimasa  Institute  of  Plasma  Physics 

Nagoya  Univ.,  Nagoya  4 64 

Takeuchi  Satoshi  Faculty  of  Engineering 

Yamanashi  Univ.,  Kofu  400 

Tanaka  Motohiko  Institute  for  Fusion  Theory 

Hiroshima  Univ.,  Hiroshima  730 
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Terashima  Yoshinosuke 

Institute  of  Plasma  Physics 
Nagoya  Univ.,  Nagoya  464 

Wakatani  Masahiro  Plasma  Physics  Laboratory 

Kyoto  Univ.,  Gokasho,  Uji,  Kyoto  611 

Yabe  Takashi  Institute  of  Laser  Engineering 

Osaka  Univ.,  Suita  565 

Observers 

U.S.A. 

Horton  Wendell  University  of  Texas,  Austin 

Institute  for  Fusion  Studies 
Austin,  Texas  78712 

Sadowski  Walter.  D.O.E. 

KOREA 

Kim  Soo  Yong  Korea  Institute  of  Technology 

Department  of  Physics 

400  Kusong-dong  Chung-gu  Taejon-shi 

Chung  chong  nam-do  Korea 

EGYPT 

Makar  Malak  N.  Assuit  National  Univ. 

Dept,  of  Mathematics 
Assuit,  Egypt 
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GENERALIZED  REDUCED  MHD,  AND  FULL  MHD  CALCULATIONS  WITH 
SEMI-IMPLICIT  TECHNIQUES 


A.  Aydemir 

University  of  Texas  Austin 
Institute  for  Fusion  Studies 
Austin,  Texas  78712 


The  Four-Field  model  of  Hazeltine,  et  al  1 2 3  ,  generalizes  the  well 

known  reduced  MHD  (RMHD)  equation  to  include  parallel  and  perpendicular 
compressibility,  diamagnetic  dift  frequencies,  and  pressure  gradients  in 
the  parallel  Ohms  law.  This  model  is  applied  to  the  study  of  drift¬ 
tearing  modes  in  tokamaks;  newly-found  Alfven-resonant  modes-  are 
discussed.  Applications  of  recently-developed  semi-implicit  techniques^ 
for  efficient  treatment  of  shear  Alfven  waves,  ion  sound  waves,  and  semi- 
col  lisional  terms  are  presented. 

semi-implicit  techniques  are  discussed  in  the  context  of  non-reduced, 
full  MHD  equations  also.  Fully  toroidal,  linear  and  nonlinear  studies  of  m 
-  1  modes  in  tokamaks  are  presented  and  discussed  in  terms  of  their  appli¬ 

cation  to  fast  saw-tooth  crashes. 

(1)  R.D.  Hazeltine.  et  al,  Phys.  Fluids,  28,  2466  (1985). 

(2)  A.V.  Aydemir.  et  al,  Phys.  Fluids,  30,  4(1987). 

(3)  D.S.  Harred,  and  D.D.  Schnack,  J.  Comp.  Phys.  (1986). 


SEMI -IMPLICIT  METHODS  FOR  3-D  MHD  COMPUTATIONS 


D.D.  Schnack 

Science  Applications  International  Corp. 
San  Diego,  Calif.  92121 


Nonlinear  MHD  systems  often  evolve  on  time  scales  that  are  long 
compared  to  those  associated  with  the  fastest  normal  modes  of  the  system. 
Stability  restrictions  placed  on  explicit  temporal  approximations  may 
result  in  uneconomically  small  time  steps,  while  the  fully  implicit  treat¬ 
ment  of  nonlinear  terms  requires  iteration,  or  results  in  unacceptably 
large  matrices.  Recently,  a  new  class  of  methods  for  solving  the  time- 
dependent  MHD  equations,  based  on  an  algorithm  developed  for  long-range 
weather  simulation,  has  been  introduced.  These  semi-implicit  methods  allow 
very  large  time  steps,  yet  avoid  the  complexity  and  large  memory  require¬ 
ments  associated  with  implicit  methods.  In  the  semi  implicit  methods  for 
MHD  new  terms  are  added  to  the  time-discretized  equations  that  do  not 
affect  the  consistency  of  the  solution,  yet  provide  a  simple  and  efficient 
means  of  enhancing  stability.  This  method  is  unconditionally  stable  with 
respect  to  all  Alfven  modes,  and  consequently  permits  such  large  time  steps 
that  accuracy  becomes  the  most  important  consideration  in  the  chioice  of 
step  size.  The  semi-implicit  method  may  be  implemented  using  a  variety  of 
different  time  advance  schemes  and  choices  for  the  semi-implicit  terms. 
The  proper  choices  can  lead  to  improved  efficiency  and  accuracy  for  long¬ 
time  scale  3-d  nonlinear  MHD  computations.  Examples  of  such  simulations  of 
fusion  devices  (tokamaks  and  RFPs)  and  the  solar  corona  are  presented. 


«■ 
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NONLINEAR  EVOLUTION  OF  FREE  BOUNDARY  MODE 
IN  A  TOKAMAK 


Gen  ichi  KURITA,  Tomonori  TAKIZUKA,  Masaf umi  AZUMI 
and  Tatsuoki  TAKEDA 

Department  of  Thermonuclear  Fusion  Research 
Naka  Fusion  Research  Establishment 
Japan  Atomic  Energy  Research  Institute 
Naka  machi ,  Naka  gun,  Ibaraki  ken 


Nonlinear  MHD  calculations  of  the  m/n-2/1  free  boundary  mode  including 
tearing  mode  in  a  cylindrical  tokamak  are  carried  out  by  taking  account  of 
the  parallel  diffusion  of  resistivity.  When  1.75<qaSl.9  for  qo/qQ=0.5  ( 
qo  is  the  safety  factor  at  the  magnetic  axis  and  qt.  is  that  at  the  plasma 
surface  ),  the  plasma  column  shrinks  with  the  elliptic  deformation,  the  value 
of  qn  is  decreased  in  time,  the  plasma  becomes  stable  against  the  m/n  2/1 
mode,  and  finally  damping  oscillation  is  observed.  The  vacuum  bubbles 

inside  the  plasma  are  formed  for  1 ,9<q„<2.  due  to  free  boundary  kink  mode 
and  for  2.0squS2.25  due  to  surface  tearing  mode,  respectively.  The 
\  saturation  of  magnetic  islands  are  observed  for  2.25<qa<2.4  due  to  surface 

tearing  mode  with  plasma  expansion  and  for  higher  qa  values  due  to  usual 
tearing  mode  with  fixed  plasma  surface  position,  respectively.  Interaction 
between  the  plasma  and  material  limiter  causes  the  shrinkage  for  all  the 
unstable  values  of  qa  to  free  boundary  mode.  When  qo  is  nearly  equal  to 
or  larger  than  unity,  the  plasma  shrinks  rapidly  and  qa  can  be  reduced  less 
than  unity  below  which  the  m/n-1/1  kink  mode  becomes  unstable.  This  plasma 
shrinkage  is  a  candidate  of  the  major  disruptions  in  the  tokamak  discharge 
with  qa  around  2. 


1.  INTRODUCTION 


In  the  tokamak  discharges  with  q. •  nearly  equal  to  2,  some  major 
disruptions  are  observed  which  limit  the  maxi mum  plasma  current  [1 ,2j. 
The  m/n  2/ 1  tearing  mode  (  m  and  n  are  the  polo: dal  and  toroidal  mod 
numbers,  respectively  )  is  considered  to  play  on  importnat  role  in  this  n:  i jor 
disruption  process.  This  disruption  process  has  been  studied  by  numerical 
calculations  [3,4,5].  The  m/n  2  1  free  boundary  kink  mode  is  also 
considered  to  play  an  importnat  role  in  the  major  disruption.  The  following 
scenario  has  been  supposed:  The  plasma  deforms  elliptically  due  to  the 

growth  of  the  kink  mode  and  the  deformation  is  saturated  by  the  negative 
surface  current.  With  the  dissipation  of  the  surface  current  due  to  th' 
plasma  limiter  interaction  [6]  or  due  to  high  electric  resistivity  near  the 
plasma  surface  [7],  however,  the  deformation  continues  to  grow  and  the  current 
disruption  is  caused.  A  new  scenaroi  of  the  disruption  has  been  presented 
by  Kurita  et  al .  [8]  by  means  of  numerical  calculaiton  of  the  mil  2  1 
free- boundary  kink  modes,  where  the  resistivity  evolution  including  parallel 
diffusion  is  considered.  The  plasma  is  deformed  by  the  interaction  with 
the  limiter.  This  shrinkage  may  cause  the  major  disruption.  The  plasma 
deformation  by  the  m  n  2  1  surface  tearing  mode  has  been  also  presented  by 
Kur  ita  et  al.  [9],  In  this  report  we  review  these  results  and  present 
detailed  results  of  numerical  calculations.  In  section  2,  basic  equations 
are  described.  The  linear  stability  of  free  boundary  modes  and  tearing 
mode  including  resistivity  equation  with  finite  parallel  diffusion  is 
analysed  in  section  3.  The  results  of  nonlinear  evolutions  are  shown  in 
section  4.  and  summary  and  discussion  are  given  in  section  5. 

2.  BASIC  EQUATIONS 

As  basic  equations,  we  employ  the  single  helicity  reduced  set  of 
resistive  MHD  equations  of  a  low  beta  cylindrical  plasma  including  the 
resistivity  evolution  equation.  Fourier  expanded  equations  in  poloidal 
and  toroidal  directions  are  written  as  follows; 


,  =  ['k,4>l*/n  +  El  Vm'/n'Jm'/i i"  E  OaO^nO 

01  »=»+*"n=n'+n“ 

—f1  =  Ch >$],»/«  n  +KlAq«/r,  , 

Uw/n  ”  A$a/n  > 

Jn/n  *  A'l'm/n  +  2 Ra/n&ttO&nO  * 


ix, Yw „  -  E  E  f- 


»'/n' 


»=»  n=n  +n 


d.\  ir"/ n" 
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(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


*„/„  is  the  helical  poloidal  magnetic  flux,  $„/„  the  stream  function,  /?„/„ 
the  resistivity,  J„/n  the  current  density,  1',/n  the  vorticity,  £7  the  electric 
field  at  the  wall,  6U  the  Kronecker  delta  and  /4/n  is  the  ratio  of  m/n,  where 
m  and  n  are  the  Fourier  mode  numbers  of  our  interest.  In  these 
equations,  the  uniform  plasma  density  is  assumed,  and  the  time  is  normalized 
by  the  poloidal  Alfven  transit  time  t pa=Bi/JpR  (Bf  is  the  toroidal  magnetic 
field,  p  the  plasma  density,  and  R  the  major  radius).  Other  normalization 
factors  are  JpR/(Btb2)  (b  is  the  shell  radius)  for  77  and  JpR/iBiR2)  for 

K|. 

The  resistivity  is  assumed  to  follow  the  same  equation  as  that  for  the 
electron  temperature.  The  parallel  diffusion  coefficient  of  resistivity. 
k,,  and  perpendicular  one,  k±,  are  assumed  to  be  uniform  for  simplicity. 
The  parallel  gradient  of  r,  is  defined  as  (v;  :•  ..  To  calculate 


the  free  boundary  problem,  we  use  the  'pseudo  vacuum’  model,  where  the  vacuum 
is  replaced  by  the  plasma  with  high  resistivity.  This  method  has  been 
successfully  applied  to  nonlinear  simulations  of  free  boundary  modes 
[10,11].  The  above  set  cf  nonline'*-  equations  is  solved  by  the 

predictor-corrector  time  integration  scheme.  The  diffusion  terms  in 
Eqs.(l)  and  (3)  are  approximated  by  the  implicit  representation.  These 
implicit  parts  of  nonlinear  calculation  spend  almost  CPU  time  of  the 
computer.  The  equation  for  resistivity,  Eq.(3),  including  the  diffusion 
term  is  solved  by  a  mapping  method.  In  this  method,  the  radial  plasma 
displacement,  S,  is  solved  instead  of  resistivity  equation,  that  is, 

^  -  [S,p]  -  Ki([T.h.'I')  +  ^(V?)2)  +  K.A *S  (7) 

The  additional  term,  (hV^) (f<i?)2.  is  appeared  by  the  variable  transformation 
of  r?-S,  and  plays  a  very  important  role  in  the  resistivity  evolution  of 
free-boundary  case.  Obtaining  new  S,/„  and  Fourier-composed  3(r,0),  we 
calculate  p(r,0)  by  a  mapping, 


r?(^,0)=p(5 t=o.0)  ;  St=o=r. 


(8) 


After  expandeing  rj(r ,0)  in  Fourier  series,  we  can  calculate  the  right  hand 
side  of  Eq.(l),  and  can  advance  the  time  step.  By  this  procedure,  the 
resistivity  is  always  kept  positive  at  all  grid  points,  which  is  necessary 
to  integrate  Eq.(l)  in  time  without  numerical  instability  [10], 

3.  LINEAR  STABILITY  ANALYSIS 


In  this  section,  we  investigate  the  effect  of  parallel  diffusion  on 
the  linear  stability  of  tearing  mode  and  free  boundary  modes;  kink  and 
surface  tearing  modes  [10,12].  From  Eqs.(l)~(6),  the  following  linearized 
reduced  set  of  resistive  MHD  equations  is  derived, 


yM  =  FM  -  ,  (9) 

-  -F<&  +  He qM  +  T)Jeq  ,  (10) 

yn  -  -  +  K,F(  )  -  on 

F  -  ^-(m-nq)  ,  (12) 


where  q  is  the  safety  factor,  subscript  ’eq*  means  equilibrium  quantities, 
and  the  time  derivatives  are  replaced  by  the  growth  rate,  y.  The 
perpendicular  diffusion  is  negligibly  small  in  general.  Since  the  singular 
surface  (  F-0  )  does  not  exist  in  the  plasma  region  for  the  kink  mode  ( 
qa<m/n  ),  y4>  st-F$  holds  and  Eq.(ll)  becomes 

(  1  +  ^  )(  yr?  +  1$  )  «  0  .  (13) 


This  relation  implies  that  the  parallel  diffusion  scarcely  affects  the  kink 
mode  in  the  linear  stage.  The  ratio  of  the  plasma  radius,  a,  to  the  wall 
radius,  b,  is  a/b=0.66.  The  resistivity,  r?eq(r),  is  inversely  proportional 
to  deq(r)  with  )?eq (0) “  1 0  6  and  ??eq(b)'=l.  The  profile  of  current  density  is 
chosen  as 
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Jeq(r)  -  {  Jeq(O)  -  Je<4(6,>  j{  1  -  (r/U)3  56  r  +  Jeq(b'),  (14) 

for  Osrga,  and  Je<!(r)  =  Je q(b)  =  Jcq(0)j?eq(0)/y?*-q(b)  for  a<rsb .  The  ratio, 

qo/q«,  is  0.5  for  this  current  profile,  and  the  value  of  JCq(0)  is 
determined  by  the  value  of  qo.  In  Fig.l,  the  graph  of  linear  growth  rate 
versus  qa  for  three  values  of  ks  is  shown  by  three  solid  curves.  The 

dependence  to  k,  is  small  in  f ree- boundary  kink  mode  region  (  qa<2)  as 
indicated  by  Eq.(13).  Numerical  calculations  of  the  eigenvalue  problem 
support  above  prediction.  The  growth  rate  of  the  kink  mode  is  smoothly 
connected  to  that  of  the  "surface  tearing  mode"  in  the  q0=2  region  [12], 
The  growth  rate  for  higher  k„  value  is  higher  than  that  of  lover  k,  at  qu 
around  2.  This  relation  is  inverted  in  the  surface  tearing  mode  region, 
and  the  growth  rate  for  higher  k,  becomes  smaller  than  that  of  lower  case. 
For  sufficiently  large  k,  value,  the  growth  rate  of  tearing  mode  with 
resistivity  perturbation  tends  to  that  without  resistivity  perturbation 
(q= const .)  shown  by  the  broken  curve  in  the  figure.  This  result  is 
consistent  to  that  in  Refs  .(13, 14],  The  relatively  large  value  of  the  linear 
growth  rate  at  qa=2.0  is  attributed  to  finite  and  relatively  large  plasma 
resistivity  near  the  singular  surface,  details  of  which  will  be  presented 
elsewhere.  The  parallel  diffusion  of  resistivity  becomes  very  important 
in  the  nonlinear  phase,  especially  in  the  phase  of  the  interaction  between 
plasma  and  limiter. 

4.  NONLINEAR  CALCULATIONS 

In  this  section,  we  carry  out  nonlinear  calculations  of  the  m/n-2/1 
free  boundary  mode  including  tearing  mode  for  the  cases  (4. a)  without  a 
limiter  and  (4.b)  with  a  limiter.  The  initial  profile  of  current  density 
is  given  by  Eq.(14)  (  qo/qJo<r0.5  ).  Total  plasma  current  is  assumed  to 
be  constant  in  time.  The  resistivity  at  10  is  determined  as 
q(r)-Eu'/J(r)  with  q(0)  =  10‘6  and  r?(b)  =  l.  We  choose,  in  the  following 
calculations,  the  initial  plasma  radius,  ao,  as  ao/b^O. 66.  The  parallel 
diffusion  coefficient,  k,  is  set  102,  and  Ki=108.  The  nonuniformity  of  q 
on  a  magnetic  surface  disappears  within  the  time  interval  of  about  0.1  due 
to  the  diffusion  of  k«=102.  This  value  of  K,  corresponds  to  the  following 
actual  parameters;  toroidal  magnetic  field  Bt=4  T,  major  radius  R=1  m, 
electron  temperature  Te=l  keV,  and  plasma  density  n=10“°  nf3.  Number  of 
Fourier  components,  M,  and  radial  meshes,  Nr,  are  typically  M=10  and 
Nr=200. 

4. a  Case  without  limiter 

The  nonlinear  evolutions  without  limiter  are  studied  at  first  for 
various  initial  values  of  qaO=qa(I=0);  (a)  qcO=l  .85,  (b)  qo0=l  .95,  (c) 
qoO-2.05,  (d)  qao=2. 15,  and  (e)  qao=2.5.  Figures  2  and  3  show  the  time 
evolutions  of  ^-contour  for  qo0=1.85  and  qao=1.95,  respectively.  The  bold 
line  in  the  figure  represents  crowded  resistivity  contours  which  correspond 
to  the  approximate  position  of  the  plasma  surface.  For  qoo=1.85,  the 
elliptic  deformation  grows  with  the  shrinkage  of  the  plasma  column  ( 
0<l£l00  ).  Since  the  plasma  current  is  constant,  this  shinkage  makes  qa 
value  smaller  from  1.85  to  1.25,  and  the  plasma  becomes  linearly  stable 
against  the  m/n^2/l  kink  mode.  Finally  the  damping  oscillation  of  the 
shrunk  plasma  is  observed  (  t>W0  ).  The  radius  of  each  shrunk  plasma  with 
circular  cross  section  is  O.95a0  for  qao-l  .75  ,  0.85ao  for  qa0=l  .8  and  0.80ao 
for  qaO^l-SS.  The  shrinkage  of  the  plasma  column  is  caused  by  the  following 
processes.  An  )?  contour  in  a  resistive  plasma  crosses  ^-contours  near  the 
plasma  surface  due  to  the  convection,  and  the  plasma  periphery  connected  with 
the  vacuum  region  is  drastically  cooled  by  the  parallel  thermal  conduction. 

It  is  to  be  noted  for  the  case  of  Ki^O  that  the  plasma  area  is  conserved 
and  the  saturation  state  with  elliptic  deformation  can  be  realized  (10). 
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On  the  other  hand,  for  larger  q^o  value  (  qao  1  .95  ),  the  vacuum  bubbles  are 
formed  by  the  free  boundary  kink  mode,  as  shown  in  Fig. 3,  even  for  such  a 
decreasing  current  profile  (  qo/Qah o=0.5  ).  The  hot  plasma  flows  out  into 
the  'pseudo  vacuum'  region  along  the  magnetic  fild  line  and  the  vacuum  region 
penetrates  into  the  plasma  to  form  the  vacuum  bubbles.  It  was  shown  by 
Rosenbluth  et  al .  that  the  saturation  state  of  the  ideal  kink  mode  is  the 
elliptic  deformation  for  the  parabolic  current  profile  [15].  The  vacuum 
bubbles  were  found  to  be  formed  by  the  surface  tearing  mode  in  a  resistive 
plasma  for  qa0>2  and  qo/qah-o=0.5  [10],  The  bubble  formation  by  the 
free  boundary  kink  mode  for  the  initial  condition  of  qa^l  .85  and  qo=1.2 
has  been  found  by  Dnestrovskii  et  al .  using  the  'heating'  model  for  the 
transport  of  r)  in  which  the  convection  of  plasma  boundary  is  not  considered 
[11).  The  transition  from  the  formation  of  vacuum  bubbles  to  the 
shrinkage  with  elliptic  deformation  occurs  between  the  q^  values  of  1.85 
and  1.9  for  this  current  profile.  The  current  profile  determines  this 
transition  point,  but  the  detailed  mechanism  or  the  criterion  of  this 
transition  is  not  clarified  yet. 

When  qao  becomes  larger  than  2,  the  formation  of  vacuum  bubbles  by  the 
surface  tearing  mode  is  observed  even  for  finite  k,  value  just  the  same  as 
for  the  case  of  no  k(  [10],  This  bubble  formation  occurs  in  the  following 
process;  First,  the  plasma  boundary  touches  the  magnetic  islands  formed 
inside  of  the  plasma  due  to  the  horizontal  plasma  flow.  Because  of  finite 
large  Ki,  the  cold  plasma  in  the  vacuum  flows  into  the  hot  plasma  along  the 
magnetic  field  line  instantaneously  and  the  magnetic  islands  with  plasma 
current  are  put  to  the  vacuum  region,  being  separated  from  the  plasma  hot 
core.  At  this  point,  the  singular  surface  is  also  put  in  the  vacuum  and 
after  that  the  vacuum  bubbles  are  formed  just  the  same  process  as  that  in 
free  boundary  kink  mode  (see  Fig. 4). 

We  observe  the  magnetic  island  saturation  with  plasma  expansion  for 
higher  qao  values  2.05SqaoS2.3  as  shown  in  Fig. 5.  The  transition  from 
formation  of  vacuum  bubbles  to  saturation  due  to  magnetic  islands  occurs 
about  qao=2.05  for  the  plasma  of  this  current  profile  with  k^IO2.  To  see 
this,  we  rewrite  Eq .(11)  by  elimination  of  by  use  of  Eq.(10)  as  follows; 

(J+k,F2]>/  -  ^(*-k,F4>)  .  (15) 

Sign  of  resistivity  perturbation  at  the  plasma  surface,  r?(a),  determines 
whether  vacuum  bubbles  are  formed  or  not,  so  the  condition  for  the  formation 
of  vacuum  bubbles  can  be  written  by  the  following  condition  for  k, ; 


$  q 

m-nq  r  a 


(16) 


In  the  case  of  Ki=0,  since  m-nq  <  0  and  ♦  >  0  at  plasma  surface  r=a,  this 
condition  becomes  to  be  $  <  0,  just  the  same  as  the  condition  of  surface 
tearing  mode  [10,12].  So,  in  this  case,  the  plasma  always  forms  the  vacuum 
bubbles  by  the  surface  tearing  mode,  which  is  consistent  to  our  former 
results  [10],  In  the  case  of  k,— this  condition  cannot  be  satisfied  and 
in  this  case  plasma  always  saturate  to  form  the  magnetic  islands.  When 
qao  value  becomes  much  larger,  and  the  unstabel  mode  changes  from  surface 
tearing  mode  to  tearing  mode,  the  usual  magnetic  islands  saturation  with 
fixed  plasma  surface  position  is  realized  in  the  simulation  (see  Fig. 6). 

4.b  Case  with  limiter 

In  this  subsection,  the  effect  of  the  limiter  on  the  nonlinear  evolution 
of  free-boundary  kink  mode  and  surface  tearing  mode  is  studied.  We  assume, 
for  simplicity,  that  the  plasma  is  surrounded  by  the  limiter  for  all  poloidal 
and  toroidal  angles.  The  limiter  radius,  bi ,  is  chosen  as  bj/ao=1.03  for 
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all  cases.  The  value  of  Ki  in  the  plasma  region  is  10°,  while  that  in 
the  vacuum  region  is  set  10  4  for  numerical  reasons.  Other  parameters  are 
the  same  as  those  for  the  case  without  the  limiter. 

At  first,  we  investigate  the  effect  of  limiter  to  the  case  of  m  r2  1 

free  boundary  kink  mode;  1 .75Sqtlo<2.0.  We  show  the  time  evolutions  of  qa 

and  the  internal  inductance,  1,,  for  the  cases  with  (  solid  line  )  and  without 
(  broken^line  )  the  limiter  in  Fig. 7  for  qao  1  .85  (a^  and  qao  1  95  (b),  where 
U  -  I  ( i)2dS  /  (v\I/)a2S  (  S  is  the  plasma  area  and  a  denotes  the  plasma 

surfaci  )  and  its  initial  value  is  1.  The  bold  line,  shown  upper  side  of 

the  figure,  denotes  the  time  duration  when  the  plasma  contacts  with 
limiter.  In  the  shrinkage  phase  with  the  constant  total  current,  the  values 
of  qo  is  unchanged  in  contrast  with  the  decrease  of  qa  value,  the  positive 
skin  current  flows  near  the  plasma  surface,  and  the  value  of  Z ,  becomes 
small.  In  the  case  of  qoO'1.85,  the  final  value  of  qa  is  about  1.3  and  that 
of  Z,  is  about  0.7,  respectively,  for  both  cases  with  and  without  the 
limiter.  After  the  plasma  surface  touches  the  limiter  (t>55)  the  plasma 
is  shrinking  as  the  ellipticity  is  increasing.  In  contrast  to  this  case, 
the  time  evolutions  of  qa  and  Z,  for  quo  1  .95  are  quite  different  between 
two  cases  with  and  without  limiter.  Figure  8  shovs  the  time  evolutions  of 
and  q  contours  for  the  case  of  qao=1.95,  in  the  region  of  formation  of 
vacuum  bubbles.  The  deformation  of  plasma  is  much  restricted  by  the  limiter 
shown  by  a  row  of  small  rectangles  in  the  figure.  After  the  plasma  touches 
the  limiter  (Zs40),  it  is  hardly  sharpened  by  the  limiter,  shrinks  rapidly 
and  finally  goes  into  the  stable  state  to  m/n  2/1  free  boundary  kink  mode. 
When  the  plasma  becomes  stable  against  the  m/n  2/1  kink  mode,  the  plasma  is 
detatched  from  the  limiter  and  and  the  damping  oscillation  begins  (fs 80), 
which  is  clearly  seen  by  the  magnetic  energy  evolution  of  each  mode  shown 
in  Fig. 9.  Since  the  separatrix  always  crosses  the  limiter,  the  vacuum 
bubbles  cannot  be  formed  and  plasma  periphery  is  cooled  as  the  same  as  the 
case  of  qaoSl.9.  As  shown  by  solid  lines  in  Fig. 7(b),  the  minimum  value 
of  qa  becomes  to  be  less  than  one  and  the  value  of  Z,  becomes  about  0.6  in  this 
case.  To  demonstrate  the  plasma  shrinkage  and  the  fluttening  of  average 
current  density,  the  radial  current  profiles  of  m/n  0/0  mode  before  (t=0) 
and  after  (t=187.5)  the  shrinkage  are  shown  in  Fig. 10.  This  phenomenon  is 
much  different  from  that  without  the  limiter;  the  formation  of  large  vacuum 
bubbles.  The  evolutions  of  qa  and  Z*  of  that  case  are  shown  by  broken  lines 
in  Fig. 7(b)  and  the  evolution  of  4'  and  rj  are  shown  in  Fig. 2,  respectively. 

Next,  we  investigate  the  limiter  effect  in  the  case  of  surface  tearing 
mode,  2.sqaoS2.45.  In  Figs.  11  and  12,  time  evolutions  of  4<  and  q  contours 
are  shown  for  qoO-2.05  and  qo0=2.15,  respectively.  In  both  cases,  when  the 
singular  surface  is  put  into  the  vacuum  region,  the  magnetic  islands  with 
ths  plasma  current  are  also  put  there.  This  helical  current  flowing  in 
the  magnetic  islands,  which  is  separated  from  the  hot  plasma  core,  makes 
the  plasma  in  the  equilibrium  state  stable  to  m/n^2/l  kink  mode  and  the 
plasma  is  detached  from  the  limter  once  (tsl30  in  Fig.  11  and  t>230  in 
Fig. 12).  The  helical  current,  however,  disappears  due  to  the  interaction 
with  the  limiter  and  the  plasma  becomes  unstable  to  m/n=2/l  kink  mode  again, 
and  the  same  process  develops  as  that  in  the  free-boundary  kink  mode.  The 
values  of  qa  become  nearly  equal  to  1  in  both  cases. 

Figure  13  is  the  stability  diagram  in  the  (qa,qo/qa)  plane  for 
approximated  current  profile  to  that  realized  in  the  nonlinear  calculations 
with  the  limiter,  where  the  hatched  region  denotes  the  unstable  one.  The 
trajectories  of  qo^const.  for  qoo=I  .85  and  qa<H  .95  are  also  depicted,  which 
correspond  to  the  results  of  nonlinear  calculations.  In  this  stability 
calculation,  the  plasma  radius  is  determined  from  the  condition  of  constant 
total  current,  that  is,  a=aov/daAfao  with  ao/b-0.66.  It  is  easily  seen  from 
the  figure  that  higher  qoO  values  for  the  same  qo/qa  value  reslut  in  lower 
qa  value  in  the  final  stable  state.  Final  q„  values  of  nonlinear 
calculations  for  both  q^  cases  are  shown  in  the  figure  by  open  circles. 
Because  of  its  inertial  effect,  the  final  q„  values  decrease  lower  than  the 
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margmaly  stable  values.  It  should  be  noted,  however,  that  the  numerical 
results  here  are  obtained  from  the  single  helicity  calculation  and  in  an 
actual  plasma  with  qa  less  than  unity,  the  most  dangerous  m  n  1,  1 
free  boundary  kink  mode  becomes  unstable,  which  can  easily  lead  the  current 
disruption.  The  existence  of  this  m/n  1,1  free  boundai  y  kink  mode  in  a 
tokamak  is  already  detected  in  T-10  experiment .and  the  relation  to  major 
disruption  is  also  suggested  in  Ref. [16],  To  show  the  whole  disruption 
process,  multiple  helicity  calculation  must  be  carried  out. 

It  is  also  easily  seen  from  Fig.  13  that,  if  the  plasma  can  pass  the 
dangerous  zone,  qa  of  about  2,  without  large  deformation  of  plasma  surface 
in  some  way,  the  m/n  2  1  free-boundary  kink  mode  becomes  not  so  dangerous 
and  the  shrinkage  due  to  the  instability  makes  the  value  of  qu  a  littiel 
bit  smaller;  just  the  intermediate  one,  1 .3S  qa  Si. 7.  where  the  plasma  is 
stable  to  both  m/n -2/1  kink  mode  and  m/n  1/1  surface  tearing  and  kink 
modes.  This  can  explain  the  reason  of  realizing  the  very  low  qu  value  in 
a  few  tokamak  experiments,  once  q£1  is  lowered  blow  2  [17]. 

4. c  Negative  surface  current 

In  this  subsection,  the  evolution  of  negative  surface  current,  which 
is  considered  to  play  an  important  role  in  the  major  disruption  process,  is 
studied.  Figure  14  shows  the  time  evolution  of  the  maximum  value  of 
negative  surface  current  near  top  or  bottom  of  plasma  poloidal  plane  in  Fig. 8 
for  qa(r- 1 . 95  without  1  imiter  (broken  line)  and  with  limiter  (solid  line'  . 
The  negative  surface  current  developed  to  suppress  the  instability  does  not 
disappear  in  spite  of  the  interaction  between  the  plasma  and  limiter,  but 
even  grows  rapidly  during  the  contact  and  disappears  suddenly  when  the  plasma 
reaches  the  stable  state  and  is  detached  from  the  limit er.  This  fact  is 
different  from  the  prediction  of  Kadomtsev(6j  or  Zakharov(7],  but  it  is 
plausible  because  the  plasma  surface  is  always  formed  just  inside  of  limiter, 
the  plasma  is  unstable  during  the  interaction  and,  furthermore,  the  negative 
surface  current  grows  almost  exponentially. 

5.  SUMMARY  AND  DISCUSSION 

In  this  section,  we  summarize  the  results  of  nonlinear  calculations. 
The  parallel  diffusion  of  resistivity  plays  a  crucial  role  in  the  evolution 
of  free-boundary  modes,  especially  in  the  case  with  the  limiter;  The  plasma 
shrinks  due  to  the  diffusion,  while  the  plasma  area  and  the  value  of  qa  is 
conserved  in  time  without  the  diffusion.  The  final  states  of  nonlinear 
m/n=2/l  free  boundary  modes  evolutions  without  the  limiter  are  classified 
into  following  four  cases  according  to  the  value  of  initial  qa,  qa0. 

1 .  Saturation  of  magnetic  islands  without  plasma  surface  deformation 

2.  Saturation  of  magnetic  islands  with  pla  la  expansion 

3.  Saturation  due  to  formation  of  vacuum  bubbles 

4.  Stable  state  of  shrunk  plasma 

The  ranges  of  qoO  are  4£qa0>2.3,  Z.3ZqM>2.05,  2.05>qa0Sl.9  and 

l-9£qoo£l.75  for  cases  1  ,  2,  3  and  4,  respectively,  which  are  values  for 
calculation  parameters  in  this  report,  and  depend  on  the  current  profile 
and/or  the  value  of  Kt .  The  boundary  between  cases  1  and  2  almost  coincide 
with  that  between  tearing  and  surface  tearing  modes,  and  that  of  cases  2 
and  3  is  given  by  Eq.(16).  As  for  the  boundary  between  cases  3  and  4, 
however ,  the  mechanism  is  not  clarified  yet,  and  we  can  determine  this 
boundary  only  by  the  nonlinear  calculations  at  present. 

When  the  plasma  is  surrounded  by  the  limiter  placed  near  the  plasma 
surface,  the  above-mentioned  final  states  of  cases  2  and  3  becomes  completely 
different  one,  that  of  case  4,  and  in  this  case,  all  unstable  free  boundary 
modes  with  plasma  surface  deformation  go  into  stable  state  of  shrunk  plasma 
to  m/n-2/1  kink  mode  .  In  some  cases,  however,  the  minimum  value  of  qa 
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becomes  even  less  than  1 ,  where  the  most  dangerous  m/ii  1/1  free  boundary  kink 
mode  is  unstable.  New  disruption  scenario  of  tokamak  discharges  with 
qa= 2  is  described  briefly  in  the  following; 

1.  When  qu~2,  m/n-2/1  free-boudary  kink  mode  becomes  unstable. 

2.  Plasma  shrinks  due  to  the  interaction  with  the  limiter. 

3.  If  qoKl.  this  shrinkage  makes  qa  value  to  be  less  than  or  nearly  equal 

to  1 . 

4.  Then  m/n=l/l  free-boundary  kink  mode  becomes  unstable,  which  causes 
the  major  disruption. 

This  process  occurs  in  the  neighborhood  of  the  broken  line. in  the  stability 
diagram  shown  in  Fig. 13.  It  can  be  concluded  that  the  condition,  qcr^l , 
is  required  to  pass  the  dangerous  zone  of  qa=2  in  tokamak  discharges, 
because,  by  the  interaction  with  the  limiter,  qa  decreases  nearly  equal  to 
1  for  qo=l  and  the  m/n=l/l  kink  mode  is  destabilized.  The  condition  of 
qo<l  for  realization  of  low  qa  (qaSl.5)  discharges  is  the  same  as  that  in 
Ref. (7],  but  the  mechanisms  are  completely  different  each  other.  The 
nonlinear  calculations  of  multiple  helicity  are  now  being  carried  out  and 
the  complete  results  of  this  disruption  scenario  will  be  presented 
elsewhere. 
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FIG.  1  Linear  growth  rate,  7.  versus  safety  factor  at  the  plasma  surface 
Qa  for  various  values  of  parallel  resistivity  diffusion  coefficients  k,  for 
current  profile;  (\-(r/a)3  56)2  Parameters  are  chosen  as  a/b=0.66 

and  j?(b)/)?(0)  =  106.  Values  of  k,  are  1 05 , 1 02  and  10°  for  three  solid  lines. 
Broken  line  is  growth  rate  without  resistivity  perturbation,  Eq.(ll). 


T I  ME  =  185-00 


TIME  =  207-50 


FIG. 2  Time  evolutions  of  ^  contours^  and  plasma  surface  (bold  line)  for 
current  profile,  J(r)-J(0) (l-Cr/a)3'56)2.  with  qao  1.85.  Plasma  shrinks  with 
elliptic  deformation. 


FIG. 3  Time  evolutions  of  ♦-contours  and  plasma  surface  for  qo0=1 .96 
Formation  of  large  vacuum  bubbles  is  observed. 
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FIG. 4  Time  evolutions  of  ^-contours  and  plasma  surface  for  qoO-2.05 
Formation  of  vacuum  bubbles  is  observed  for  plasma  of  qas2  with  finite  K,. 


FIG. 7  Time  evolutions  of  qa  and  internal  inductance  lt  for  (a)  q^l  .85 
and  (b)  <ja0=1.96.  Results  of  calculation  with  and  without  limiter  are  shown 
by  solid  and  broken  curves,  respectively.  Bold  line,  shown  upper  side  of 
figure,  denotes  the  time  duration  when  plasma  contacts  with  limiter. 


FIG. 8  Time  evolutions  of  ♦-contour  and  plasma  surface  for  qoO'1.95. 
Limiter  is  placed  at  r  l  ,03cio.  The  plasma  touches  limiter  at  38  and  is 

detached  from  it  at  t= 79. 
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FIG. 9  Tim-.’  ev : I « •  L : u.nr.  of  magnetic  energy  of  each  mod*, 
limiter.  Bold  lin-*.  slv-wn  upper  side  of  f r*  .  duiv>Le? 
when  the  plasma  c  .'‘lit  acts  with  linn  ter. 
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FIG  10  Radial  current  profiles  of  m/n  0/0  mode  before  (t-0;  broken  line) 
and  after  (t-187.5;  solid  line)  plasma  shrinkage  for  qao=1.95  Plasma 
radius,  a/b0.66  initially,  becomes  a/1) -0.485  through  this  shrinkage. 


TIME  =  155-00  TIME  =  167. SO 


FIG. 11  Time  evolutions  of  ^-contour  and  plasma  surface  for  surface  tearing 
mode  with  qao=2.05.  After  helical  current  disappears  due  to  interaction  with 
limiter,  the  same  phenomenon  occurs  as  that  in  Fig. 8. 
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FIG. 12  Time  evolutions  of  ^-contour  and  plasma  surface  for  surface  tearing 
mode  with  qc0=2.15.  Almost  the  same  phenomenon  occurs  as  that  in  Fig. 11. 
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FIG. 13  Stability  diagram  of  m,/n=2/l  free  boundary  kink  mode  for 
approximated  current  density  profile  to  that  realized  in  nonlinear 
calculations.  Hatched  region  denotes  unstable  one,  and  two  broken  lines 
of  trajectory  correspond  to  nonlinear  calculations  with  limiter  for 
qo0=l  .85  and  qao=1.95.  Black  and  open  squares  denote  start  and  final 
positions,  respectively.  Trajectory  of  qo0=~2  (qo=l)  is  shown  by  dotted 
line.  Plasma  radius  is  determined  from  condition  of  constant  total  current, 
i.e.,  a=ao/\Ala/qaO  with  ao/b=0.66. 

Current  density  profiles  are  determined  from  following  formula, 

j (r)  =  jo(  (1  - (r/ao)3  56  )  2+dexp(-2( (r~a+p)/V)2)  ]{ 1  -(r/a) v}2. 

Height  of  skin  current,  A,  is  determined  from  condition  of  constant  total 
current  for  fixed  position  of  skin  current,  p=0.1ao,  and  width  of  that, 
y=0.05ao.  Value  of  v,  which  denotes  current  gradient  near  plasma  surface, 
is  chosen  to  be  several  ten. 
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FIG. 14  Time  evolutions  of  maximum  value  of  negative  surface  current  near 
top  or  bottom  of  plasma  poloidal  plane  in  Fig. 8.  This  is  the  same 
calculation  as  in  Fig. 7(b).  Bold  line,  shown  upper  side  of  figure,  denote 
the  time  duration  when  the  plasma  contacts  with  limiter.  Broken  line 
denotes  result  without  limiter. 
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A  high  precision  MHD  code,  which  has  the  fourth-order  accuracy  for  both 
the  spatial  and  time  steps,  is  developed,  and  is  applied  to  the  simulation  studies  of 
two  dimensional  driven  reconnection.  It  is  confirm  that  the  numerical  dissipation 
of  this  new  scheme  is  much  less  than  that  of  two-step  Lax-Wendroff  scheme.  The 
effect  of  the  plasma  compressibility  on  the  reconnection  dynamics  is  investigated  by 
means  of  this  high  precision  code. 
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§1.  Introduction 


The  magnetohydrodynamic  (MHD)  simulations  have  given  a  great  contribu¬ 
tion  to  the  researches  of  the  nuclear  fusion  and  astrophysical  plasmas.  The  widely 
used  numerical  schemes  of  the  MHD  simulation  codes  are  categorized  to  two  types 
of  schemes,  i.e.  Fourier  and  finite  difference  schemes.  The  former  is  useful  for  the 
problems  that  have  a  periodic  boundary  condition,  and  the  latter  is  favorable  for  the 
free  boudary  problems,  the  typical  example  being  space  plasma  phenomena.  Two- 
step  Lax-Wendroff  scheme^  is  one  of  the  most  widely  used  finite  difference 
schemes.  This  scheme  time-centres  the  integration  by  defining  temporary  or  inter¬ 
mediate  values  of  the  dependent  variables  at  the  half  time  steps,  and  this  is  a  explicit 
method  that  has  the  second-order  accuracy  in  the  spatial  and  time  steps.  Because  of 
the  stability  and  the  convenience  for  the  implementation,  two-step  Lax-Wendroff 
scheme  is  much  powful  for  analyzing  the  global  MHD  dynamics.  In  this  scheme, 
however,  large  numerical  diffusion  occurs  in  short  wavelength  components,  and 
hence  this  method  is  not  so  appropriate  for  the  problem  in  which  short  wavelength 
components  are  important,  for  example,  turbulence  and  shock  problems. 

In  this  research,  we  develop  a  high  precision  MHD  code  that  has  the  fourth- 
order  accuracy  in  both  the  spatial  and  time  steps,  and  apply  that  to  the  simulation  of 
two  dimensional  driven  reconnections.  It  is  widely  known  that  there  are  two  types 
of  reconnections,  one  being  tearing  mode  reconnection  and  the  another  being  driven 
reconnection.  Sato  and  Hayashi  in  detail  investigated  the  dynamics  of  driven  recon¬ 
nection  by  means  of  two-step  Lax-Wendroff  simulation  code,  and  they  found  that 
driven  reconnection  leads  to  accelaration  of  plasma  as  high  as  the  local  Alfven 
speed,  and  that  slow  MHD  shocks  is  formed  just  downstream  of  the  separatrix. 
They  also  revealed  that  once  reconnection  proceeds,  its  ultimate  fate  no  longer 
seems  to  be  dependent  on  the  resistivity,  but  is  largely  controlled  by  the  boundary 
conditions. 
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The  first  objective  of  this  paper  is  to  reconfirm  the  above  conclusions  by 
means  of  high  precision  code.  The  second  one  is  to  show  the  effectiveness  of  the 
high  precision  code  by  quantitatively  comparing  the  simulation  results  by  two  dif¬ 
ferent  codes.  And  the  another  objective  is  to  investigate  the  influence  of  the  plas¬ 
ma  compressibility  on  the  reconnection  dynamics  by  using  of  the  high  precision 
code. 


§2.  Numerical  Scheme  and  Simulation  Model 

We  employ  the  fourth-order  finite  difference  approximation  to  the  spatial 
derivative 

(f£)j  =  j  =  3,  4,...,  JV,_2 

where  the  simulation  region,  x=  -Lx/2~Lx/2,  is  implemented  on  a  set  of  Nx  mesh 
points  xf,  Xj=  -Ly2  +  0'-l)Ax,  &x-Lx/(Nx- 1).  For  the  other  spatial  dimensions, 
the  same  fashion  is  employed.  And  the  time  integral  is  carried  out  by  means  of  the 
fourth-order  Runge-Kutta-Gill  method.  The  time  integral  scheme  for  eq. 
3f/dt=F(f)  is  given  as  follows, 

Af°  =  F(f„°)-A/ 
f.1  =  f„°+c1-Af° 

Af1  =  F^-A/ 

f*2  =  f„,  +  c2-(Af1-Af°) 

Af,1  =  G2-Afn0+S2-Af! 

Af2  =  F(f*2)-Af 

f,3  =  f„2  +  c3-(Af2-Afl) 
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Af.2  =  <23  Af„'  +  S3  Af2 

Af3  =  F(f*3)Af 
f„+i°  =  f„3+c4-(Af3-2-Af.2) 

where  cj=  1/2,  c2=  1—  1/V2,  c3  =  1+  1A/2,  c4=  1/6,  ()2  =  l~3c2,  £)3=l-3c3, 

S2=2c2,  S3=2c3,  and  n  denotes  the  time  step. 

The  boundary  condition,  the  initial  condition  and  the  basic  equations  are  just 
same  as  Sato  and  Hayashi’s  ones.  The  initial  distributions  of  the  magnetic  field 
Bx(z),  the  current  j(z),  and  the  pressure  p(z),  are  chosen  to  be 

Bx(z)  =  flx0tank(z),  j(z)  =  7osech2(*)>  (i) 

p(z)  =  p0sech2(z)  +  P0,  (2) 

where  Po~Bx02/ 2,  and  P0  is  constant  for  space.  (But  P0  is  usually  zero,  if  not  indi¬ 
cated  differently.)  The  simulation  region  is  a  square  box  surrounded  by  four 
planes,  two  being  parallel  to  the  x  axis  and  placed  at  z  =  ±2  (normalized),  and  the 
other  two  being  parallel  to  the  z  axis  and  placed  at  x-  ±3.  We  inject  the  plasmas 
with  a  given  mass  flux,  magnetic  energy  flux,  and  total  energy  flux  symmetrically 
from  the  boundaries  placed  at  z  =  ±2  (input  boundary).  The  boundaries  at  x=  ±3 
(output  boundary)  are  assumed  to  be  free  boundary.  We  presume  that  the  resistivi¬ 
ty  t)  takes  the  form 

i)(J)  =  a(j~jc)2  f°r  j^jc 

0  otherwise. 

And  the  basic  equations  to  be  solved  are  the  compressible,  isotropic,  one-fluid  mag- 
netohydrodynamic  equations.  The  real  calculation  is  carried  out  by  using  of  the 
normalized  set  of  these  equations  as  same  as  in  Ref. 2. 
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§3.  Results 

Figure  1  is  the  contour  plot  of  the  magnetic  flux  at  three  different  times, 
where  Aq  is  the  maximum  speed  of  the  injected  flow  from  the  input  boundaries. 
As  same  as  Fig.l  in  Ref. 2,  the  X-type  neutral  point  and  the  separatrix  are  formed 
before  t=  12  Alfven  times.  Afterwards,  reconnection  proceeds  and  plasma  is  ac¬ 
celerated  as  high  as  local  Alfven  speed,  the  slow  shock  being  generated.  Figure  2 
shows  the  cross  section  distributions  of  j,  Bx  and  p  in  the  results  by  two-step  Lax- 
Wendroff  code  and  high  precision  code.  All  parameters  in  both  calculations  by  two 
codes  are  same.  Note  that  the  global  structures  of  these  two  simulations  are  almost 
same.  This  global  agreement  strongly  supports  the  correctness  of  the  main  result  in 
Ref.2,  i.e.  the  slow  shocks  are  formed  just  downstream  of  the  separatrix.  Further¬ 
more,  we  can  see  that  the  high  prcision  code  more  sharply  reproduce  the  shock 
structure  on  about  2-0.5.  In  the  up-stream  side  of  shock  front,  however,  the 
current  distribution  by  the  high  pricision  code  has  a  wavy  structure.  These  small 
differences  in  the  local  structure  between  two  simulations  are  caused  by  the  reduc¬ 
tion  of  the  numerical  diffusion  in  the  high  prcision  code. 

Figures  3,  4,  and  5  show  the  temporal  evolutions  of  the  electric  field  by  both 

two-step  Lax-Wendroff  and  the  high  precision  codes,  when  a,jc,  and  Aq  are 

changed.  No  remarkable  change  between  two  different  calculations  confirms  Sato 

and  Hayashi’s  conclusion  that  the  crucial  parameter  for  the  electric  field  on  the 

reconnection  point,  i.e.  for  the  reconnection  rate,  is  not  the  resitivity,  but  the  plasma 

2 

flow  on  the  input  boundary,  Aq.  The  interesting  difference  between  two-step 
Lax-Wendroff  code  and  the  high  precision  code  is  that  the  pulsation  of  the  electric 
field  is  more  clearly  observed  in  the  high  precision  code  than  in  two-step  Lax- 
Wendroff  code.  The  difference  clearly  appears  in  the  curves  of  Aq=0.2  in  Fig. 5. 
This  pulsation  is  caused  by  the  effect  of  magnetosonic  wave,  which  is  bouncing 
between  two  input  boundaries.  The  disappearance  or  tne  reduction  of  magnetoson- 
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ic  wave  pulsations  in  two-step  Lax-Wendroff  code  seems  to  be  caused  by  large  nu¬ 
merical  diffusion  in  two-step  Lax-wendroff  code.  These  results  show  that  the  nu¬ 
merical  diffusion  of  the  high  precision  code  is  less  than  that  in  two-step  Lax- 
Wendroff  code. 

Lastly,  we  investigate  the  effect  of  the  plasma  compressivility  on  the  recon¬ 
nection  dynamics  by  this  high  precision  code.  Though  driven  reconnection  have 
been  studied  numerically  by  many  authors,  the  most  of  them  employed  the  in- 
compressible  approximation.  Recently,  however,  it  is  revealed  in  the  reversed- 
field  pinch  study  that  the  rate  of  driven  reconnection  is  reduced  by  the  effect  of  the 
finite  3  of  plasma.  The  reason  is  that  the  converging  flow  into  the  X-point  can 
lead  to  the  peaking  of  the  plasma  pressure  on  the  X-point  by  the  influence  of  the 
adiabatic  compression,  and  the  peaked  pressure  on  the  X-point  acts  to  suppress  the 
converging  flow  and  to  delay  the  reconnection  process.  From  this  results,  we  can 
suspect  that  as  the  system  is  approaching  to  the  incompressible  system  the  driven 
reconnection  is  being  reduced,  since  the  incompressible  approximation  is 
corresponding  to  the  large  pressure  limit  or  large  magnetic  field  limit.  In  order  to 
confirm  this  suspection,  we  execute  some  calculations  of  several  P0  in  eq.(2).  The 
results  are  shown  in  Fig. 6,  where  the  up-  and  down-side  panels  show  the  evolutions 
of  the  plasma  density  and  electric  field  on  the  reconnection  (X-)point,  respectively. 
In  this  simulation,  the  plasma  resistivity  is  settled  to  a  constant  value  (ii  =  0.01),  and 
the  injected  plasma  flow  is  adjusted  so  that  the  electric  field  on  the  input  boundaries 
is  constant  (Eq  =  0.2).  We  can  see  that  as  P0  is  increasing  the  change  of  the  density 
on  the  reconnection  point  is  reduced.  If  we  remember  that  the  uniform  density 
profile  is  unchanged  everywhere  in  the  incompressible  approximation,  we  can  con¬ 
sider  that  the  system  is  approaching  to  the  incompressible  system  as  P0  is  increasing. 
And  we  can  see  also  that  the  electric  field,  that  is  proportional  to  the  reconnection 
rate,  is  saturated  in  lower  level  for  larger  P0.  Therefore,  we  can  conclude  that  the 
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incompressible  effect  leads  to  the  reduction  of  the  rate  of  driven  reconnection. 
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2-step  Lax-Wendroff  high  precision  code 


Fig. 2.  Cross  sectional  distributions  of  the  current  (jy) ,  x 
component  of  the  magnetic  field  ( Bx ),  pressure  (/>),  and  den¬ 
sity  (p)  on  a  plane  at  x  =  2.4.  Note  the  opposite  changes  of  Bx 
and  p  (and  p )  across  the  current  layer,  thus  showing  evidence 
of  slow  shocks. 


T 


2-step  Lax-Wendroff 


Fig. 3.  Developments  of  neutral  point  electric  fields  for  two 
different  threshold  currents  jc.  Note  that  there  exists  no 
marked  difference. 
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2-step  Lax-Wendroff 


Fig.4.  Developments  of  neutral  point  electric  fields  for  three 
different  resistivity  coefficients  a.  Note  that  there  exists  no 
marked  difference. 
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Fig. 6.  Developments  of  plasma  density  and  electric  field  on  X-point  for 
four  different  Pq. 
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ABSTRACT 

Interaction  between  the  solar  wind  and  the  plasmas  produced  from 
cometary  gases  has  been  modeled  by  using  a  three-dimensional  time-dependent 
magnetohydrodynamic  (MHD)  simulation  with  the  cometary  mass  loading.  The 
model  reproduced  several  features  observed  by  the  recent  missions  to  comet 
Halley.  A  weak  bow  shock  was  located  at  3.2x1  C>5  km  in  front  of  the  comet. 
The  magnetic  field  increased  hy  a  factor  of  3-7  a ctosb  this  weak  bow  shock 
and  continued  to  increase  up  to  the  contact  surface  (BmnY/BTMTi,=6->-8)  ■  The 
plasma  temperature  increased  across  the  bow  shock  and  decreased  nearer  to  the 
comet.  IMF  lines  were  hung  up  on  the  comet  and  formed  a  long  plasma  tail  in 
which  the  lobe  field  was  quite  strong  and  B]_0t>e/BD^p-w5.  A  cold  dense  plasma 
sheet  formed  in  the  tail  and  this  thin  plasma  sheet  was  oriented  normal  to 
the  IMF  direction. 

1.  Introduction 

Interaction  between  the  solar  wind  and  a  comet  is  somewhat  different 
from  that  between  the  solar  wind  and  unmagnetized  planets  like  Venus.  The 
difference  may  occur  because  the  comet  has  a  very  small  mass  or  a  small 
gravitational  force.  Therefore,  a  large  amount  of  gas  is  evaporated  from  the 
cometary  icy  nucleus  and  extends  into  the  solar  wind  where  it  is  ionized  by 
the  photoionization  process  due  to  the  solar  ultra-violet  ray  and  the  charge 
exchange  process.  Thus  a  quite  broad  dayside  interaction  region,  which  is 
composed  of  three  discontinuities,  the  outer  shock,  contact  surface  and  inner 
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shock,  can  be  formed  around  the  comet  [Brant  and  Mendis,  1979;  Mendis  and  Ip, 
1977;  Mendis  and  Houpis,  1982;  Niedner,  1984]*  The  cometary  plasma  tail 
frequently  has  a  complicated  structure  with  narrow  rays,  kinks  and 
condensations.  Moreover  it  sometimes  shows  a  mysterious  behaviour  which  is 
known  as  a  disconnection  event  [Niedner  and  Brandt,  1978,  1979;  Niedner  et 
al.,  1981].  That  is,  the  cometary  tail  is  disconnected  and  propagates  down 
the  tail  from  the  comet. 

In  the  international  spacecraft  missions  to  comet  Halley,  the  Suisei 
probe  got  to  a  closest  approach  of  1 . 51x1 O^km  at  13:06  UT  on  6  March  1986 
[Itoh  and  Hirao,  1986],  The  plasma  flow  vectors,  the  proton  density  and  the 
estimated  directions  of  draped  magnetic  field  in  the  sun  side  strong 
interaction  region  were  demonstrated  during  the  Suisei  encounter  with  comet 
Halley  [Mukai  et  al.,  1986a,  b].  The  transition  corresponding  to  the  bow 
shock  was  clearly  shown  and  the  shape  of  the  bow  shock  was  given  by  a 
parabola.  The  shock  distance  from  the  comet  was  estimated  to  be  3*5x1C)5km  on 
the  sun-comet  line  at  the  time  of  the  closest  approach  of  the  Suisei  probe. 
The  Geotto  magnetometer  experiment  [Neubauer  et  al.,  1986]  demonstrated  that 
the  magnetic  field  was  strongly  enhanced  (B=50%60nT)  near  the  comet  and 
suddenly  disappears  in  the  cometary  center  region  or  probably  inside  the 
contact  surface.  The  many  features  of  the  strong  interaction  region  have 
been  made  clear  from  the  satellite  observations.  On  the  other  hand,  a 
disturbance  of  the  plasma  tail  of  comet  Halley  was  also  observed  by  handreds 
of  photographs  and  its  dynamics  was  discussed  in  connection  with  the 
disconnection  events  of  cometary  plasma  tail  [Saito  et  al.,  1986]. 

Schmidt  and  Wegmann  [1982]  solved  the  full  set  of  MHD  equations  to  model 
the  interaction  between  the  solar  wind  and  the  cometary  ionosphere  in  the 
region  exterior  to  a  fixed  contact  surface  by  introducing  the  cometary  gas 
production.  Two  surfaces  of  discontinuity,  the  shock  front  and  the  contact 
surface  are  reproduced  and  the  embedded  magnetic  field  induces  an  asymmetry 
in  the  plasma  flow  and  in  the  density  distribution  of  cometary  tail.  Fedder 
et.  al.  [1986  a,b]  also  simulated  the  interaction  processes  to  demonstrate 
the  overall  structure  of  the  cometary  tail.  A  spheroidal  head  and  a  long 
ribbon-like  tail  were  formed  for  a  steady  solar  wind  and  interplanetary 
magnetic  field  (IMF).  Furthermore,  predictions  for  comet  Giacobini-Zinner 
were  presented  by  the  MHD  simulation  of  comets.  The  interaction  of  the  solar 
wind  with  the  outflowing  plasma  from  a  comet  was  simulated  by  using  a  2 
dimensional  MHD  model  to  reproduce  the  contact  surface  and  the  cometary 
magnetotail  [Ogino  et  al.,  1986b].  For  a  constant  IMF,  tail  magnetic 
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reconnection  begins  to  occur  at  several  points  in  a  thin  plasma  sheet.  A 
plasmoid  with  high  plasma  density  appears  in  the  cometary  tail  and  propagates 
tailward  with  a  local  Alfven  speed.  When  the  IMF  orientation  is  reversed, 
dayside  magnetic  reconnection  occurs  at  the  subsolar  point,  and  a  larger 
disturbance  occurs  in  the  whole  tail  region  and  propagates  down  the  tail. 

In  the  present  paper,  the  interaction  processes  between  the  solar  wind 
and  comet  Halley  have  been  successively  studied  by  uBing  a  3  dimensional 
time-dependent  MHD  model  with  a  cometary  plasma  production.  The  model 
reproduced  several  features  of  the  comet  and  solar  wind  interaction  predicted 
by  earlier  theories  and  observed  from  the  recent  missions  to  comet  Halley. 

2.  Simulation Model 


The  present  purpose  is  to  study  the  interaction  between  the  solar  wind 
and  the  cometary  plasmas  by  using  a  3  dimensional  MHD  simulation  [Ogino, 
1986a]  with  plasma  production.  Therefore,  the  plasma  production  rate  is 
added  to  the  MHD  equations  as  it  was  given  by  Schmidt  and  Wegmann  [1982]. 
The  normalized  MHD  equations  are  written  as  follows, 
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stands  for  the  cometary  plasma  production,  J=VxB  is  the  current  density,  P 
the  plasma  density,  v  the  flow  velocity,  p  the  plasma  pressure,  B  the 
magnetic  field,  i=uVzv  the  viscosity,  7=5/3  the  ratio  of  specific  heats, 
7?=7?o ( T/T0 ) ~3/2  the  resistivity  and  T/T0  the  normalized  temperature.  For  the 
sake  of  convenience  the  units  of  distance ,  velocity  and  time  are  respectively 
defined  by  the  earth's  radius  (Re=6.37x10^m) ,  the  Alfven  velocity  (va=6.80x 
lO^m/s)  and  the  Alfven  transit  time  (ts=Re/vA=0.937s) .  The  other  numerical 
parameters  are  »o=0.01,  tt/Pgy=D=Dp=0.005  where  Psw  is  the  solar  wind  density, 
and  the  small  diffusion  terms  in  (la),  (lb)  and  (1c)  are  simply  added  in 
order  to  suppress  numerical  MHD  fluctuations.  For  the  present  comet 
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simulation,  vc  amd  pc  for  the  cometary  plasma  are  neglected  because  |vc|  « 
|v|  and  pc  «  p.  The  magnetic  Raynold's  number,  which  is  the  magnetic 
diffusion  time  devided  by  the  Alfven  transit  time,  S  is  larger  than  400. 

The  coordinate  system  used  in  the  3  dimensional  MHD  simulation  is  shown 
in  Figure  1,  where  x0=y0=z0=7. 7x10 5km  and  x-|=-1 . 54x1  O^km.  The  solar  wind 
accompanying  a  uniform  IMF  in  the  z-direction  flows  into  the  simulation  box 
through  the  boundary  at  x=xQ  and  begins  to  interact  with  the  cometary  plasma. 
The  mirror  boundary  is  used  at  z=0  and  y=0,  on  the  other  hand  the  free 
boundary,  where  the  spatial  derivative  of  the  physical  quantities  is  zero,  is 
used  at  z=zD,  y=y0  and  x=x-j . 

The  typical  parameters  of  the  solar  wind  are  vsw=500km/s,  nsw=15/cc, 
Taw=2x1o5°K  and  Bjj,jp=Bz=6nT  for  a  closest  approach  of  the  Suisei  probe  to 
comet  Halley  [Itoh  and  Hirao,  1986;  Mukai  et  al.,  1986a,  b].  Moreover,  the 
parameters  of  the  plasma  production  for  comet  Halley  Eire  as  follows:  the 
plasma  production  rate,  Q0=1  .OxIO^Os-"' ,  the  ionization  rate,  cf=3.3x10~^s~^ , 
the  radial  flow  velocity,  vr=1km/s,  the  ionization  distance,  Ac=vr/a 
=3. 03x1 05km  and  the  effective  cometary  mass  ratio  of  the  water  group  to  the 
proton,  mc=l6  [Mukai  et  al.,  1986a,  b].  The  MHD  equations  are  solved  by  the 
two  step  Lax-Wendroff  method  as  an  initial  value  problem  on  the  grid  points 
of  (Nx,  Ny,  Nz)=(90,  30,  30)  ,  {150,  50,  50)  or  (120,  60,  60)  except  for  the 
boundary.  The  mesh  size  is  AX=Ay=Az=4Re,  3Re  or  2Re  and  the  time  step,  At  is 
selected  as  8Ax  (=30s  for  Ax=4Re)  in  order  to  assure  the  numerical 
stability. 

3.  Simulation  Results 


Interaction  processes  between  the  solar  wind  and  comet  Halley  were 
simulated  by  using  the  observed  parameters  for  the  closest  approach  of  the 
Suisei  probe  [Mukai  et  al.,  1986a,  b]  as  was  previously  mentioned.  Stating 
from  the  initial  conditions  where  a  uniform  solar  wind  only  exists,  the  code 
was  run  until  a  quasi-steady  state  configuration  resulted.  This  required 
about  1024  time  steps  or  8.5  hours  in  real  time. 

In  Figure  2  is  shown  a  quasi-steady  state  configuration  of  the 
interaction  between  the  solar  wind  and  comet  Halley  in  the  x-z  and  x-y  planes 
for  no  uniform  IMF  where  Bjj^p=Bz=0nT,  vsw=500km/s  and  nsw=15/cc.  The  plasma 
density,  P  and  the  plasma  pressure,  p  are  depicted  by  the  contours  and  the 
flow  velocity  by  the  arrows.  A  bow  shock  is  formed  at  3*2x1o5km  as  is  shown 
in  the  plasma  pressure  and  the  velocity  flow  pattern,  and  a  cometary  dense 
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plasma  extends  toward  the  downstream.  The  plasma  flow  changes  its  direction 
toward  the  outside  across  the  bow  shock,  however  it  is  rather  parallel  to  the 
sun-comet  line  in  the  further  tail.  The  configuration  of  interaction  is 
symmetric  in  the  x-z  and  x-y  planes  because  of  no  uniform  IMF.  When  a 
uniform  IMF  of  Bu,ip=B2=6nT  is  introduced  in  the  z  direction,  the  quasi-steady 
state  configuration  is  modified  as  is  shown  in  Figure  3*  On  the  top  panel 
the  magnetic  field  is  also  depicted  by  the  arrows.  The  position  of  the  bow 
shock  moves  little  and  an  asymmetry  appears  in  the  x-z  and  x-y  planes.  The 
plasma  density  and  pressure  are  thinner  in  the  z  direction  parallel  to  the 
uniform  IMF,  which  comes  from  the  draping  effect  of  hanged  IMF  lines.  The 
magnetic  field  has  a  strongly  disturbed  feature  in  the  neibourhood  of  the  bow 
shock  and  is  again  well  ordered  in  the  tail  region.  The  field  magnitude  is 
quite  large  even  in  the  tail  lobe  (B=30nT).  It  should  be  noted  that  a  high 
density  lump  appears  at  x=-8.8x10^km  in  the  tail  and  propagates  tailward. 
The  tailward  velocity  is  about  50km/s  which  is  comparable  with  a  local  Alfven 
velocity,  v^=40km/ s . 

In  Figure  4  are  shown  the  profiles  of  the  magnetic  field,  B, 
temperature,  T,  x-component  of  velocity,  vx  and  density,  P  on  the  sun-comet 
line  where  M  is  an  effective  mass  ratio  to  the  proton.  The  weak  bow  shock  is 
formed  at  x=  3.2x105km  and  the  Mach  number  of  the  fast  magnetosonic  wave  is 
about  M=2.0  on  the  sun-comet  line,  which  is  quite  reduced  from  the  solar  wind 
value,  M=6.45.  This  is  because  that  the  plasma  temperature  increases  whereas 
the  velocity  decreases  in  the  upstream  region  of  the  bow  shock  due  to  the 
cometary  plasma  production  processes.  In  fact,  the  simulated  temperature, 
velocity  and  density  are  T/M=1 .3x10^°K,  v=350km/s  and  p=Mn=22/cc  in  front  of 
the  bow  shock.  The  temperature  increases  across  the  weak  bow  shock  and  then 
decreases  nearer  to  the  comet.  The  magnetic  field  increases  by  a  factor  of 
3.7  across  this  weak  shock  and  continued  to  increase  up  to  the  contact 
surface  (Bmax/Bji<[p^6) .  A  cold  and  dense  plasma  extends  toward  the  tail  and 
the  high  density  lump  is  seen  at  x=-8.8x1  C^km.  The  magnetic  field  also 
exhanced  on  the  comet  side  of  the  high  density  lump  at  x=-7.6x1()5km. 

In  Figure  5  are  shown  the  profiles  of  the  magnetic  field,  B, 
temperature,  T,  density,  P,  velocity,  v  and  the  direction  of  velocity,  0V  at 
x=1.7x1()5km  in  the  dayside  interaction  region.  The  sharp  gradients  at 
z=y-5x105km  for  T/M,  B,  v  and  P  correspond  to  the  weak  bow  shock.  The 
velocity  somewhat  decreases  even  in  the  outside  of  the  weak  shock  from  the 
solar  wind  value,  vsw=500km/s  due  to  the  plasma  production.  The  flow 
velocity  decreases  from  460km/s  to  260km/s  across  the  weak  shock  and 
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gradually  approaches  a  minimum  value  of  80km/s  on  the  sun-comet  line.  The 
plasma  density  and  the  magnetic  field  increase  about  2.7  times  across  the 
shock  and  continue  to  increase  toward  the  sun-comet  line.  The  flow  angle,  8V 
deviates  about  25°  from  the  sun-comet  direction  (8^180°)  inside  the  bow 
shock.  These  features  are  very  similar  to  the  Suisei  probe  observation  when 
they  are  compared  with  Figure  2  in  the  article  of  Mukai  et  al.  (1986a).  The 
plasma  temperature  increases  about  6.8  times  across  the  weak  bow  shock  and 
decreases  toward  the  sun-comet  line.  It  is  noted  that  the  magnetic  field 
strongly  fluctuates  along  the  z-direction  as  well  as  the  y-direction  inside 
the  bow  shock. 

In  Figure  6  are  shown  the  profiles  of  the  magnetic  field,  B, 
temperature,  T,  density,  P  and  velocity,  v  at  x— 8.1x1 o5km  in  the  cometary 
tail.  The  magnetic  field  is  enhanced  up  to  about  5  times  the  IMF  value  in 
the  tail  lobe  and  is  a  little  enhanced  up  to  10nT  along  the  y-axis 
perpendicular  to  the  uniform  IMF.  The  magnetic  field  does  not  extremely 
decreases  near  the  sun-comet  line  at  x=-8.1x1C)5km,  which  means  the  intrusion 
of  the  magnetic  field  into  the  cold  and  dense  cometary  plasma  tail,  however 
it  becomes  very  small  in  further  tail  for  xc-I.OxIO^xm  (see  Figure  4).  The 
temperature  and  the  velocity  decrease  at  the  center  of  tail  and  the  plasma 
density  increases  there.  The  high  plasma  density  in  the  tail  is  thin  in  the 
z-direction  or  in  the  uniform  IMF  direction  and  expands  in  the  y-direction, 
which  is  consitent  with  the  draping  picture  of  hanged  IMF  lines. 

In  Figure  7  are  shown  the  cross  sectional  patterns  of  the  plasma 
pressure,  p,  density,  p  and  the  x-component  of  magnetic  field,  Bx  in  the  tail 
for  Bz=6nT  and  OnT.  Tail  lobes  and  a  thin  cold  plasma  sheet  are  deary 
formed  by  draping  effect  of  IMF  lines  (see  Bx  component).  The  cold  plasma 
sheet  is  thin  in  the  parallel  direction  to  the  uniform  IMF,  and  the  tendency 
becomes  remarkable  in  the  distant  tail.  On  the  other  hand,  the  plasma  sheet 
is  naturally  symmetric  for  no  uniform  IMF,  Bz=0nT  as  are  shown  in  the  lower 
two  panels. 

In  Figure  8  are  shown  the  three  dimensional  configuration  of  the 
magnetic  field  lines  and  their  projection  onto  the  three  planes  for  the 
interaction  between  the  solar  wind  and  comet  Halley,  where  considerably 
draped  IMF  lines  are  only  depicted  and  the  colour  means  a  measure  of  the 
field  draping.  It  is  clearly  seen  on  the  figures  that  the  IMF  lines  hang  up 
to  the  cometary  plasma  as  was  suggested  by  Alfven  [1957]  and  slip  away  toward 
the  tail.  The  slipped  magnetic  field  lines  gradually  come  close  toward  the 
sun-comet  meridian  (plane  at  y=-0)  in  the  further  tail. 
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4»  Conclusion 

The  interaction  between  the  solar  wind  and  the  cometary  plasmas  has  been 
simulated  by  using  a  three  dimensional  MHD  model  in  order  to  compare  with  the 
recent  sattelite  observations  to  comet  Halley.  As  the  result,  several 
features  of  the  interaction  which  are  consistent  with  the  Suisei  probe  and 
the  Giotto  observations  have  been  reproduced  when  the  parameters  proposed 
from  the  satellite  observations  for  the  solar  wind,  the  IMF  and  the  plasma 
production  of  comet  Halley  were  used. 

The  simulated  position  of  the  weak  bow  shock,  and  the  spatial  variation 
of  the  plasma  density,  velocity  and  temperature  in  the  sun  side  interaction 
region  are  consistent  with  those  of  the  observational  results  from  the  Suisei 
probe.  The  dense  and  cold  cometary  plasma  sheet  extends  toward  the 
downstream  of  the  solar  wind  and  becomes  thinner  down  the  tail  in  the 
parellel  direction  to  the  upstream  IMF.  The  magnetic  field  is  enhanced  up  to 
about  3  times  the  IMF  value  across  the  weak  bow  shock  and  continues  to 
increase  up  to  6-^8  times  the  IMF  value  near  the  comet.  The  field  magnitude 
of  the  draped  IMF  is  not  small  even  in  the  tail  lobes  and  keeps  about  5  times 
the  IMF  value. 

The  MHD  fluctuations  were  enhanced  in  the  neibourhood  of  the  weak  bow 
shock.  The  plasma  production,  which  in  turn  captures  the  IMF  lines  into  the 
cometary  interaction  region,  as  well  as  the  Kelvin-Helmholts  instability  due 
to  the  velocity  shear  must  have  the  responsiblity  to  the  excitation  of  the 
MHD  waves.  A  further  simulation  is  needed  to  study  the  excitation  mechanism 
of  MHD  waves  in  the  neibourhood  of  bow  shock  and  also  cometary  tail  dynamics 
such  as  the  magnetic  reconnection  and  disconnection  events. 

Acknowledgements 

We  would  like  to  acknowledge  helpful  discussions  with  Drs.  D.A. 
Mendis,  T.  Mukai  and  T.  Terasawa.  This  work  was  supported  by  a  Grant-in 
Aid  for  Science  Research  from  the  Ministry  of  Education,  Science  and  Cultrure 
and  by  the  NASA  Solar  Terrestrial  Theory  Program  Grant  NAGW-78.  The 
simulations  were  performed  at  the  Computer  Center  of  Nagoya  University. 

References 

Alfven,  H.,  On  the  theory  of  comet  tails,  Tellus,  £,  92,  1957.  Brandt,  J.C., 


-53- 


and  D.A.  Mendis,  The  interaction  of  the  solar  wind  with  comets,  in 
Solar  System  Plasma  Physics,  Vol.  II,  edited  by  C.  F.  Kennel  et  al., 
p.  253,  North  Holland  Publ.  Co.,  Amsterdam,  1979. 

Fedder,  J.A.,  S.H.  Brecht,  and  J.G.  Lyon,  MHD  simulation  of  a  comet 
magnetosphere,  Icarus,  submitted,  1983. 

Fedder,  J.A.,  J.G.  Lyon,  and  J.L.  Giuliani,  Jr.,  Numerical  simulations  of 
comets:  Prediction  for  Comet  Giacobinni-  Zinner,  EOS,  Jan.  1 4 »  1986. 

Itoh,  T.,  and  K.  Hirao,  The  Sakigake  and  Suisei  encounters  with  comet 
Halley,  Geophys.  Res.  Lett.,  13,  817,  1986. 

Mendis,  D.A.,  and  W.-H.  Ip,  The  ionosphere  and  plasma  tails  of  comets, 
Space  Scl.  Rev.,  20,  145,  1977. 

Mendis,  D.A.,  and  H.L.F.  Houpis,  The  cometary  atmosphere  and  its  interaction 
with  the  solar  wind,  Rev.  Geophys.,  20,  885,  1982. 

Mukai,  T.,  W.  Miyake,  T.  Terasawa,  M.  Kitayama,  and  K.  Hirao,  Plasma 

observation  by  Suisei  of  solar  wind  interaction  with  comet  Halley, 
Nature.  321,  299,  1986a. 

Mukai,  T.,  W.  Miyake,  T.  Terasawa,  M.  Kitayama,  and  K.  Hirao,  Ion 
dynamics  and  distribution  around  comet  Halley:  Suisei  observation, 
Geophys.  Res.  Lett.,  13,  829,  1986b. 

Neubauer,  F.M. ,  K.H.  Glassmeier,  M.  Pohl,  J.  Raeder,  M.H.  Acuna,  L.  F. 
Burlaga ,  N.F.  Ness,  G.  Musmann,  F.  Mariani,  M.K.  Wallis,  E. 
Ungstrup,  and  H.U.  Schmidt,  First  results  from  the  Geotto  magnetometer 
experiment  at  comet  Halley,  Nature ,  321 ,  352,  1986. 

Niedner,  M.B.,  and  J.C.  Brandt,  Interplanetary  Gas  XXIII.  Plasma  tail 

disconnection  events  in  comets:  Evidence  for  magnetic  field  line 
reconnection  at  interplanetary  sector  boundaries,  Astrophys .  J . ,  223, 

655,  1978. 

Niedner,  M.B.,  and  J.C.  Brandt,  Interplanetary  Gas.  XXIV.  Are  cometary 
plasma  tail  disconnections  caused  by  sector  boundary  crossings  or  by 
encounters  with  high  speed  streams?,  Astrophys.  J.,  234,  723,  1979* 

Niedner,  M.B.,  J.A.  Ionson,  and  J.C.  Brandt,  Interplanetary  Gas.  XXVI.  On 
the  reconnection  of  magnetic  fields  in  cometary  ionospheres  at  inter¬ 
planetary  sector  boundary  crossings,  Astrophys.  J..  245,  1159,  1981. 

Niedner,  M.B.,  Magnetic  reconnection  in  comets,  Geophys.  Mono.,  30,  79, 

1984- 

Ogino,  T.,  A  Three  dimensional  MHD  simulation  of  the  interaction  of  the  solar 
wind  with  the  earth's  magnetosphere:  The  generation  of  field-aligned 
currents,  J.  Geophys.  Res.,  91 ,  6791,  1986a. 


Ogino,  T.,  R.J.  Walker,  and  M.  Ashour-Abdalla,  An  MHD  simulation  of  the 
interaction  of  the  solar  wind  with  the  outflowing  plasma  from  a  comet, 
Geophys,  Res.  Lett.,  13,  929 »  1986. 

Saito,  T.,  K.  Yumoto,  K.Hirao,  K.  Saito,  T.  Nakagawa,  and  E.J.  Smith,  A 
disturbance  of  the  ion  tail  of  comet  Halley  and  the  heliospheric 
structure  as  observed  by  Sakigake,  Geophys.  Res.  Lett.,  13,  821,  1986. 

Schmidt,  H.U.,  and  R.  Wegmann,  Plasma  flow  and  magnetic  fields  in  comets,  in 
Comets:  Gases,  Ices,  Grains  and  Plasma,  edited  by  L.L.  Wilkening,  p. 
538,  University  of  Arizona  Press,  Tucson,  Arizona,  1982. 


z 


Figure  1 .  Coordinate  system  of  the  solar  wind  and  comet  Halley  interaction 
in  the  3  dimensional  MHD  simulation. 


Figure  2.  Configurations  of  the 
solar  wind  and  comet  Halley 
interaction  in  the  x-z  and  x-y 
planes  where  Bjj^OnT, 

vsw^500km/s  and  nsw=15/cc. 
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Figure  3.  Configurations  of  the 
solar  wind  and  comet  Halley 
interaction  in  the  x-z  and  x-y 
planes  where  BjMp=6nT , 

vsw=500km/s  and  nsw=15/cc. 


Figure  4.  Profiles  of  the  magnetic 
field,  B,  temperature,  T, 
x-component  of  velocity,  vx  and 
density,  p  on  the  sun-comet 
Halley  line  where  M  is  an 
effective  mass  ratio  to  the 
proton. 
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Figure  8.  Three  dimensional  configuration  and  projection  of  the  magnetic 
field  lines  in  the  interaction  regions  of  the  solar  wind  and  comet 
Halley,  where  the  colour  means  a  measure  of  the  field  draping. 


NOVA-2:  A  Kinetic-MHD  Stability  Code 
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To  study  the  effect  of  energetic  particles  on  MHD  instabilities,  a 
kinetic-MHD  stability  code  (NOVA-2)  has  been  developed.  A  low-frequency 
drift  kinetic  model  is  employed  to  describe  energetic  particle  dynamics  and 
the  MHD  fluid  model  for  the  bulk  plasma.  NOVA-2  is  and  extension  of  the 
previously  developed  nonvariational  ideal  MHD  stability  code(NOVA)  for 
computing  stability  of  low  -n  modes  in  toroidal  geometry.  Some  numerical 
examples  of  the  energetic  particle  driven  toroidicity- induced  shear  Alfven 
waves  anf  fishbone  modes  will  be  presented. 
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A  HYBRID  VLASOV-FLUID  MODEL  WITH  KINETIC  IONS  AND  MASSLESS 

FLUID  ELECTRONS* 


John  M.  Dawson 

University  of  California.  Los  Angeles 
Department  of  Physics, 

405  Hilgard  Los  Angeles,  Calif.  90024 


We  have  developed  two-and  three-dimensional  models  with  kinetic  ions 
and  massless  fluid  electrons.  The  use  of  fluid  electrons  allows  us  to 
eliminate  the  electric  field  so  that  it  never  has  to  be  solved  for.  We  are 
then  able  to  derive  a  Vlasov  equation  for  the  particle  ions  and  obtain  the 
linear  dispersion  relation  from  it.  We  compare  the  results  from  this 
theory  with  simulation  results  from  our  model;  excellent  agreement  is 
found.  The  model  exhibits  Alfven  waves,  ion  cyclotron  waves,  low  frequency 
whistler  waves  and  ion  Bernstein  waves.  We  have  used  the  model  to  simulate 
instabilities  such  as  might  be  caused  by  the  injection  of  an  energetic 
neutral  beam  into  the  plasma. 

We  have  alec  implemented  a  two  time  step  correction  technique  into 
code  which  improves  energy  conservation  by  a  factor  of  10(3)  and  allows  us 
to  use  time  steps  which  are  an  order  of  magnitude  larger. 

*Work  done  in  collaboration  with  F.  Kazeminejad,  J.N.  Leboeuf,  and  R. 


Sydora . 
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GYROKINETIC  PARTICLE  SIMULATION  OF  FINITE-BETA  PLASMAS 


W.W.  Lee  and  T.S.  Hahm  and  R.D.  Sydora* 

Plasma  Physics  Laboratory 
Princeton  University 
P.0.  Box  451,  Princeton,  N.J.  08544 


Gyrokinetic  particle  simulation  for  finite-beta  plasmas  will  be 
discussed.  The  four  main  topics  are:  1)  the  basic  formulation  of  the 
nonlinear  gyrokinetic  equations  including  magnetic  perturbations,  2)  num¬ 
erical  properties  of  the  simulation  plasma  such  as  time  step,  grid  spacing 
and  noise  level,  3)  special  techniques  for  the  evaluation  of  the  genera¬ 
lized  Ohm's  law,  and  4)  the  relationship  between  the  gyrokinetic  equations 
and  the  MHD  equations. 

*  University  of  California,  Los  Angeles  Department  of  Physics  405 
Hilgard  Ave.  Los  Angeles,  Calif.  90024 


MACROSCALE  PARTICLE  SIMULATION  OF  KINETIC  ALFVEN  WAVES 
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1  Institute  for  Fusion  Theory,  Hiroshima  University,  Hiroshima  730,  Japan 
2  AT&T  Bell  Laboratories,  Murrray  Hill,  New  Jersey  07974,  U.S.A. 

Abstract 

Simulation  of  the  kinetic  Alfven  wave,  which  has  many  applications 
in  space  and  fusion  plasmas,  has  been  performed  by  using  macroscale 
(magnetohydrodynamic  scale)  particle  simulation  code.  Dispersion 
properties  and  wave-particle  interactions  of  the  kinetic  Alfven  wave 
(such  as  Landau  damping  and  electron  acceleration)  are  successfully 
obtained  in  the  simulation. 

1.  Introduction 

The  kinetic  Alfven  wave  is  the  Alfven  wave  for  which  wave-particle 
interactions  are  important1-3.  This  wave  has  received  much  attention 
recently  in  connection  with  particle  acceleration  along  the  auroral  field 
lines4,5  and  plasma  heating  in  fusion  plasmas6.  The  kinetic  Alfven 
wave  can  also  be  an  active  agent  to  heat  the  plasma  in  the  solar  corona7 
and  in  Jupiter.  Moreover,  the  structure  of  the  auroral  arcs  seem  to  be 
determined  partly  by  this  wave8  since  the  latitudinal  scale  of  the  arcs 
(—lOOifcm)  is  comparable  to  the  ion  Larmor  radius. 

The  existence  of  the  kinetic  Alfven  wave  has  been  suggested  only 
for  limited  experimental  observations  in  space  and  in  the  laboratory. 
The  Pc  5  magnetic  pulsations  observed  on  geostationary  satellites9  show 
small  azimuthal  wavelength  300-  lOOOifcm  and  an  azimuthal  (west¬ 
ward)  propagation  speed  of  the  order  of  lOkm/sec.  The  wavelength 
and  azimuthal  phase  speed  of  the  kinetic  Alfven  wave  (k./kx«l)  are 
theoretically  scaled  by  2-irp,  and  (k./kx)vA,  respectively,  where  the  typi¬ 
cal  proton  Larmor  radius  p,  on  the  geostationary  orbit  is  lOOifcm  and  the 
Alfven  speed  vA~1000km/sec.  In  a  small  tokamak  experiment  at 
Lausanne10,  they  observed  low  frequency  density  fluctuations  and  spa- 
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tial  damping  of  the  wave  which  could  be  the  kinetic  Alfven  wave. 

Numerical  simulation  of  Alfven  waves  were  restricted  to  a  scale 
length  comparable  to  the  Debye  length  \e  when  the  conventional  parti¬ 
cle  code  was  employed.  This  is  because  the  parallel  electric  field  of  the 
wave,  which  appears  as  a  consequence  of  the  finite  ion  Larmor  radius, 
must  be  solved  properly;  but  the  space  and  time  scales  of  the  kinetic 
Alfven  waves  are  in  the  magnetohydrodynamic  scales.  So  far,  simula¬ 
tion  of  plasma  heating  by  the  antenna-launched  Alfven  wave  was  made 
using  a  long  time  scale  but  Debye  length  scale  Darwin  particle  code1 1 . 


2.  Simulation  model  and  common  parameters 

Two  types  of  simulations  of  the  kinetic  Alfven  wave12  are  presented 
here  by  using  a  newly  developed  macroscale  particle  simulation  code13 
which  enables  us  to  follow  individual  particle  dynamics  in  the  MHD 
scales,  namely,  in  the  scale  length  much  larger  than  the  Debye  length. 
In  this  code,  low  frequency  electromagnetic  electric  and  magnetic  fields 
are  solved  by  eliminating  high  frequency  oscillations  such  as  the  light 
modes;  the  scalar  potential  electric  field  is  solved  by  eliminating  Lang¬ 
muir  oscillations.  This  is  achieved  by  using  the  decentered  time  dif¬ 
ferential  scheme  which  can  be  coded  in  semi-implicit  algorithm13  or  in 
full-implicit  algorithm.  Both  ions  and  electrons  are  treated  as  particles 
with  finite  mass  whose  motions  are  governed  by  the  relativistic  version 
of  the  equations  of  motion.  The  electron  motion  along  the  magnetic 
field  line  is,  therefore,  properly  treated  in  this  code. 

The  simulation  system  is  assumed  doubly  periodic  in  the  x  and  z 
directions.  An  ambient  magnetic  field  is  applied  in  the  z  direction. 
The  velocity  is  normalized  by  the  speed  of  light  c,  the  time  by  to',1, 
and  the  electric  and  magnetic  fields  normalized  as  eE/mecb)pe  and 
eB/mec(>ipe,  respectively.  The  number  of  the  grids  is  32x64  with  the 
system  length  Lx  =  50c/(ope  and  L:  =  400c/a>pe  (except  Lx=l00c/upe  for 
the  run  in  the  next  section  and  run  A  in  Table  1).  This  system  size 
corresponds  to  nearly  250\ex2000\e.  32768  particles  are  used  for 
each  ion  and  electron  species  with  the  particle  size  of  a  =  3c/upe.  The 
parameters  used  for  the  background  thermal  plasma  are  the  ion  beta 
3  =  0.36,  the  ratio  of  the  electron  cyclotron  frequency  to  the  plasma  fre¬ 
quency  <ace/tope=l  (i.e.,eBo/mec<ape=  1),  the  ion  to  electron  temperature 
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ratio  T,/Te-4,  and  the  mass  ratio  ml/me  =  50.  Then  the  Alfven  speed 
becomes  vA  =  0.14c  and  the  ion  Larmor  radius  p1  =  v,/a)c,~3 c/o)pe  where 
v,  =  1/2  is  the  ion  thermal  speed. 


3.  Simulation  results:  Kinetic  Alfven  wave  eigenmode  in  a  thermally 

near  equilibrium  plasma 

Let  us  first  describe  the  simulation  which  successfully  detects  the 
kinetic  Alfven  wave  eigenmode  in  a  thermally  near-equilibrium  plasma 
by  means  of  power  spectrum  analysis.  As  the  initial  condition,  the 
plasma  particles  are  loaded  as  quietly  as  possible  without  giving  any 
drifts  and  density  perturbations.  Figure  1  shows  the  measured  power 
spectrum  for  the  y  component  of  the  magnetic  field  as  functions  of  the 
perpendicular  wavenumber  kx  and  the  frequency  co  for  fixed 
L.  =  2TT/(400c/ti>/,J.  The  amplitude  of  the  spectral  power  is  plotted  in  a 
logarithmic  scale  above  (to  the  right  of)  each  baseline  for  descrete  kx 
values.  Here  for  the  spectral  analysis,  the  maximum  entropy  method 
has  been  emplyed.  The  peaks  in  Figure  1  show  the  presence  of  the 
eigenmode  in  the  low  frequency  range,  i.e., a>/<oa<0.1.  The  measured 
frequency  increases  with  the  perpendicular  wavenumber  when  kx p,  be¬ 
comes  comparable  to  unity.  The  other  small  peaks  at  the  higher  fre¬ 
quencies  a)~a)c,  and  2<DCi  are  noises  due  to  ion  gyromotion.  The  line  in 
Figure  1  corresponds  to  the  theoretically  derived  wave  frequency  for 
the  kinetic  Alfven  wave.  The  dependence  of  this  measured  eigenmode 
frequency  on  the  perpendicular  wavenumber  shows  good  agreement 
with  the  theory. 


4.  Simulation  results:  Finite  amplitude  kinetic  Alfven  wave  and 
associated  particle  acceleration 

As  the  second  type  simulation,  the  propagation  of  the  finite  ampli¬ 
tude  kinetic  Alfven  wave  and  associated  nonlinear  effects  are  present¬ 
ed.  Four  runs  are  made  by  systematically  changing  the  perpendicular 
wavenumber  kx  with  the  parallel  wavenumber  k.  being  fixed  (Table  1). 
To  set  up  a  monochromatic  plane  wave  initially,  the  density  and  current 
perturbations  are  calculated  using  an  electromagnetic  kinetic  dispersion 
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solver.  The  perturbation  assumes  the  form  of  8n~sin»J/,  Sv^sinijj, 
8vv~cos\|»,  and  8v;~sin4>  where  ty  =  kxx  +  k:z.  The  amplitudes  of  these 
initial  perturbations  are  given  in  Table  1.  To  represent  the  density  per¬ 
turbation  which  is  characteristic  of  the  kinetic  Alfven  wave,  the  particle 
spacing  is  modified  in  the  particle  loading.  The  drift  velocities  of  the 
aforementioned  forms  are  given  to  the  particles  in  superposition  to  the 
random  thermal  spread. 

The  scalar  potential  field  4>,  the  x  and  z  components  of  the  electric 
field  (EX,E.)  and  the  magnetic  field  By  are  shown,  respectively,  in  Fig¬ 
ure  2(a),  (b)  and  (c)  for  t/rA  =  1.8  of  run  B  where  ta  =  2i T/k.vA  is  the 
Alfven  wave  period.  It  is  seen  that  the  loaded  wave  propagates  in  the 
direction  almost  perpendicular  to  the  ambient  magnetic  field  keeping 
well  the  initial  wave  pattern  in  the  two-dimensional  magnetic  plasma. 
The  major  fields  are  the  x -component  of  the  electric  field  Ex  and  the  y- 
component  of  the  magnetic  field  By.  A  small  but  finite  parallel  electric 
field,  a  unique  property  of  the  kinetic  Alfven  wave,  exists  as  seen  in 
Figure  2(b)  owing  to  the  finite  Larmor  radius  effect  (this  will  be  dis¬ 
cussed  below). 

The  y-component  of  the  magnetic  field  By  sliced  at  constant  x  posi¬ 
tions  for  runs  A,  B  and  C  are  shown  against  the  positions  z  and  time  t 
in  the  bird’s-eye  view  plots  of  Figure  3(a),  (b)  and  (c),  respectively. 
For  run  A  with  *vp,  =  0.19,  the  loaded  wave  propagates  with  the  parallel 
phase  velocity  close  to  the  Alfven  speed  (see  also  Table  1).  As  the  per¬ 
pendicular  wavenumber  kx  increases  for  run  C,  the  parallel  phase  velo¬ 
city  becomes  appreciably  large,  as  seen  in  Figure  3(c). 

The  phase  velocity  of  the  waves  in  the  direction  parallel  to  the  am¬ 
bient  magnetic  field,  ta/k.,  is  plotted  in  Figure  4  against  the  perpendicu¬ 
lar  wavenumber  which  is  normalized  with  the  ion  Larmor  radius,  i^p,. 
As  kx pi  increases,  the  parallel  phase  velocity  increases  parabolically  as 
shown  by  solid  circles  in  the  figure.  The  solid  line  shows  the  theoreti¬ 
cal  dispersion  relation  for  the  kinetic  Alfven  wave. 

Shown  with  the  open  circles  in  Figure  4  are  the  damping  increments 
of  the  loaded  waves.  The  numerical  damping  rate  of  the  wave  is  es¬ 
timated  to  be  yv/k.vA<3x  10~3  at  most,  which  is  negligibly  small  ex¬ 
cept  for  run  A.  For  run  D  with  kxpt~\.\,  the  damping  increment  is  so 
large  that  the  loaded  wave  decays  substantially  in  a  few  wave  periods 
[(2'n-y/ti))o/>(~0.5].  The  dashed  line  in  the  figure  shows  the  theoretical 
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damping  increment  of  the  kinetic  Alfven  wave.  Both  the  measured 
parallel  phase  velocity  and  the  damping  increment  of  the  loaded  waves 
are  in  excellent  agreement  with  the  theory. 

The  damping  of  the  kinetic  Alfven  wave  occurs  due  to  the  Landau 
damping.  Under  the  parameters  used  in  this  study,  the  damping  is  pri¬ 
marily  due  to  the  electrons.  The  initial  Maxwellian  distribution  of  the 
electrons  parallel  to  the  ambient  magnetic  field  in  Figure  5(a)  evolves 
into  a  non-symmetric  one  with  respect  to  v. -0,  which  has  a  plateau  on 
the  positive  side  around  v.~t o/Jfc-(~vf)  as  shown  in  Figure  5(b).  How¬ 
ever,  the  locations  of  both  edges  of  the  velocity  distribution  remained 
unchanged.  The  parallel  electron  effective  temperature  continues  to  in¬ 
crease,  especially  when  the  wave  amplitude  is  large.  On  the  other 
hand,  ion  parallel  temperature  stays  almost  the  same.  The  electron 
heating  rate  is  found  to  be  proportional  to  the  damping  increment  of 
the  wave  for  runs  A  to  D,  which  is  scaled  as  <37',;/dr~16-yefi  where  eB  = 
8B?/16tt  is  the  magnetic  field  energy. 

The  electrons  are  monotonically  accelerated  towards  the  direction  of 
the  wave  propagation.  The  acceleration  is  not  appreciable  for  run  A. 
For  runs  B,  C  and  D,  the  acceleration  occurs  in  the  early  time  when  the 
electric  field  energy  ef  =  &E~/16tt  is  large,  i.e.,eE/Te> 2xl0-3.  At  the 
end  of  run  D,  the  total  increase  in  the  electron  parallel  velocity 
amounts  to  nearly  10%  of  the  initial  electron  thermal  speed.  The  in- 
creasein  the  parallel  momentum  density  of  electrons  is  found,  with  an 
accuracy  of  10%,  to  be  in  agreement  with  the  decrease  in  the  total  wave 
momentum  density,  (8B'/8tt)(L/uj).  This  fact  supports  the  validity  of 
nonlinear  results  obtained  by  using  the  macroscale  particle  simulation 
code. 

The  electron  acceleration  and  heating  are  attributed  to  the  small  but 
non-vanishing  parallel  electric  field.  In  order  to  estimate  the  torque  on 
the  electrons,  the  parallel  electric  field  is  averaged  with  the  electron 
density  as  a  weight,  i.e.  ,<(- e)E.bne>/<ne> .  (The  normal  average 
field  <£->  vanishes  exactly  due  to  the  periodicity.)  The  weighted 
average  electric  field  becomes  ~lxlO“5e£0  for  run  C  when  the  elec¬ 
tron  acceleration  is  occurring  and  is  consistent  with  the  time  rate  of  in¬ 
crease  in  the  parallel  electron  momentum  (note  the  initial  wave  ampli¬ 
tude  Exk(t= 0)~3x  10~2fl0).  This  results  from  the  occurrence  of  a  posi¬ 
tive  correlation  between  the  electron  motion  and  the  electric  field,  i.e., 
the  phase  difference  between  the  electron  density  perturbation  and  the 
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electric  field  becomes  within  90  degrees.  Using  the  measured  values, 

<(  — e)E-&ne>/<nf>~ lx  10-5eB0  and  (<£->)' C~  1x10 “4fl0,  the 

phase  angle  difference  is  estimated  to  be  78  degrees.  (Note  that 
<sim}/sin(i)i  +  8)>  =  cos8/2  where  the  average  is  made  over  «|».)  Similar¬ 
ly,  the  energy  transfer  rate  to  the  electrons  from  the  wave  is  estimated 
by  measuring  <(  —  e)neve:E:>/<ne> .  This  quantity  becomes  positive 
definite,  i.e.,  4x10 ~6ceB0  for  run  C,  only  during  the  electron  heating. 
These  facts  quantitatively  verify  the  electron  acceleration  and  heating 
observed  in  the  simulation. 


5.  Conclusion 

A  successful  application  of  the  newly  developed  macroscale  particle 
simulation  code  was  shown  for  the  kinetic  Alfven  wave,  which  has 
many  applications  in  space  and  astrophysical  plasmas.  The  detailed 
wave  characteristics,  the  dependences  of  the  frequency  and  the  Landau 
damping  on  the  perpendicular  wavenumber,  were  examined  and  excel¬ 
lent  agreements  were  found  between  the  simulation  and  theoretical 
prediction.  Some  fundamental  nonlinear  interactions  of  the  kinetic 
Alfven  wave  with  the  particles  (parallel  acceleration  of  the  electrons) 
were  also  found. 

The  present  simulation  work  assumed  a  periodic  (closed)  system. 
In  the  space  environment,  however,  particle  acceleration  (heating)  oc¬ 
curs  in  an  open  geometry.  Moreover,  the  Alfven  speed  changes  its 
magnitude  in  an  inhomogeneous  plasma,  for  example,  along  the  mag¬ 
netic  field  line  of  the  earth.  In  that  case,  particles  resonating  with  the 
wave  may  be  accelerated  and  convected  toward  the  direction  of  the 
wave  propagation.  Consequently,  the  distribution  function  of  the  parti¬ 
cles  may  evolve  into  a  bi-Maxwellian  one  rather  than  formation  of  the 
plateau.  These  specific  problems  of  the  kinetic  Alfven  wave  will  be  re¬ 
ported  in  the  future. 
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Table  1 

The  initial  perturbation  of  density,  bn,  and  the  drift  velocities,  Sv^,8\J,8vJ,  for 
the  ions  (electrons)  in  the  left  (right)  column  (upper  part).  The  measured  phase 
velocity  (frequency)  (to /k:)  and  the  damping  increment  -y  of  the  wave  (lower  part). 
kx  and  kt  are,  respectively,  the  perpendicular  and  parallel  wavenumber;  p,  is  the  ion 
gyroradius;  v,  is  the  Alfven  speed. 


Run  A 

Run  B 

Run  C 

Run  D 

6  n'n0 

(0.069  .  0.0-1) 

(0.126.  0.129) 

(0.219,  0.227) 

(0  "5,  0.288) 

S.yc 

(  2.5xl0-3.-5.4xl0-4) 

(  24x  10_3,~6.4x  10-4! 

(  2.3x  10"3,  -6.'x  10~4) 

(  2.2x  10_3,--6.'x  10~4) 

6iy  c 

10-2,-:.8xl0-2) 

(-2.4x  10~2,~2.8x  10"2) 

(-  1.6x  10-2,-3.lx  10-2) 

(-8.2x  10~3,-3.6x  10-2) 

8i -Jc 

(  =0,  1.2X10"2) 

(  =0.  2.4X10'2) 

(  =0,  4.7X10-2) 

(  =0,  6.9X10-2.) 

tan  76.0 

82.9 

86.4 

87.6 

0.19 

0.38 

0.76 

1.13 

1.03 

1.07 

1.41 

1.98 

rk:vA 

(7.2=3.6)xl0~3 

(1.3  =  0.7)xl0-2 

(5.9  =  0.9)x  10~2 

(1.6=0. 2JX10-1 

► 

► 

r 


A 
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Figure  captions 


Fig.  1  The  power  spectrum  for  the  y -component  of  the  magnetic 
field  in  a  thermally  near-equilibrium  plasma.  The  abscissa  is  the 
perpendicular  wavenumber  normalized  with  the  ion  Larmor  ra¬ 
dius,  and  the  ordinate  is  the  frequency  ni/<ope.  The  power  ampli¬ 
tudes  are  plotted  in  a  logarithmic  scale  above  (to  the  right  of)  each 
baseline  with  the  104  range  between  two  adjacent  baselines. 


Fig.  2  The  scalar  potential  field  <f>,  the  x  and  z -components  of  the 
electric  field  ( EX,E. ),  and  the  y -component  of  the  magnetic  field 
By  for  t/xA  =  \.%  of  run  B.  The  maximum  amplitudes  are  (a) 
e$/mec2~  0.22,  (b)  eEx/mecu>pe~'i.5x.  10~2,  and  (c) 

eBy/mec(0p =0.17. 


Fig.  3  The  y-component  of  the  magnetic  field  sliced  at  constant  x 
positions  plotted  against  position  z  and  time  t  for  runs  A,  B  and  C 
from  top  to  bottom,  respectively. 


Fig.  4  The  parallel  phase  velocity  (solid  circles)  and  the  damping 
increment  (open  circles)  of  the  initially  loaded  waves  of  runs  A  to 
D.  The  solid  and  dashed  lines,  respectively,  show  the  theoretical 
parallel  phase  velocity  and  the  damping  increment. 


Fig.  5  The  electron  velocity  distribution  function  parallel  to  the 
ambient  magnetic  field  at  (a)  f=0  and  (b)  t=1.8TA  of  run  C.  The 
scale  in  the  ordinates  is  the  same  in  (a)  and  (b). 
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PARTICLE-PARTICLE  PARTICLE-MESH  COTE  FDR  NONIDEAL  HIGH  DENSITY 
PLASMA  AND  ITS  APPLICATION  ON  RAYLEICH-TAYLOR  INSTABILITY 

IN  ICF  PLASMA* 


Katsunobu  Nishihara 

Institute  of  Laser  Engineering,  Osaka  University 
Suita,  Osaka  565 


In  an  inertial  confinement  fusion  plasma,  the  Coulomb  coupling  con¬ 
stant,  T  =  (Ze)~  ctkgT  ,  becomes  of  the  order  of  one,  where  a  is  the  Wigner- 
Seitz  radius  defined  as  a~  (3  4zni  '■ 1/3  .  In  such  a  strongly  coupled  plasma 
(T  ^  1 )  ,  the  Coulomb  interaction  can  not  be  treated  as  perturbations  and 
the  system  begins  to  exhibit  features  qualitatively  different  from  those  in 
an  ideal  collisionless  plasma.  To  study  dynamics  of  a  strongly  couples 
plasma,  we  have  developed  three-dimensional  two-component  particle-particle 
particle-mesh  code.  The  code  calculates  exact  Coulomb  forces  among  parti¬ 
cles  within  a  short  distance,  while  the  PIC  method  is  used  for  long 
distance  forces.  We  have  observed  the  pair  distribution  functions  and  the 
result  has  been  found  to  agree  quite  well  with  the  hypernetted  chain  theory 
and  Thomas-Fermi  model . 

We  use  the  code  to  study  effect  of  a  contact  potential  on  the 
Rayleigh-Taylor  instability  between  two  plasmas  with  different  charge 
state.  The  preliminary  results  will  be  discussed. 

*Work  done  in  collaboration  with  H.  Furukawa  and  M.  Kawaguchi 
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TRIPIC:  Triangular-Mesh  Particle-in-Cell  Code  for  LIB  Diode 


Shigeo  Kawata ,  Masami  Mat sumo to  and  Yukio  Masubuchi 

The  Technological  University  of  Nagaoka 
Nagaoka,  Niigata  940-21 ,  Japan 


A  particle-in-cell  code  using  triangular  meshes  (TRIPIC  code)  is 
described  to  simulate  the  motion  of  charged  particle  in  an 
electromagnetic  field  self-ccnsistenly  in  some  complicated  computational 
space  region.  TRIPIC  is  fitted  for  the  simulation  of  LIB  (Light  Ion 
Beam)  diode,  for  example. 


1 .  Introduction 

Up  to  now  many  particle  simulations  have  been  dene  [  1  -5  ]  and  many 
particle-in-cell  (PIC)  codes  [2,6,7]  have  been  developed  to  study  a 
self-consistent  interaction  between  charged  particles  and 
electromagnetic  field.  In  many  codes  usually  rather  regular  space 
meshes  are  employed  to  describe  the  computational  regions.  On  the  other 
hand,  recently,  we  have  a  problem  of,  for  example,  intense  LIB  diode 
simulation  in  which  the  computational  space  region  has  a  nonuniform  or 
curved  boundary.  TRIPIC  code  employs  triangular  space  meshes  to 
simulate  such  the  problem.  The  paper  presents  the  description  for 
TRIPIC  code. 

In  order  to  solve  the  electromagnetic  field  on  triangular  meshes, 
we  have  a  finite  element  method  and  finite  differencial  one  which  is, 
for  example,  developed  by  Winslow  [8].  Here  we  focus  on  the  latter 
method.  Winslow  developed  the  numerical  method  in  nonuniform  triangular 
meshes  to  solve  the  Poisson  type  of  equation.  SUPERFISH  code  [9]  was 
also  developed  by  the  similar  method.  A  descendant  of  these  two  codes 
is  TRIDIF  code  [10]  which  is  a  time -dependent  code  to  solve  a  magnetic 
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diffusion.  Based  on  this  method  the  Poisson  equation  and  the  equation 
of  magnetic  vector  potential  are  solved  in  TRIPIC  code. 

In  TRIPIC  code  a  particle  pusher  is  solved  by  the  Buneman  scheme. 
The  interaction  between  particles  and  space  meshes  are  accomplished  by  a 
simple  weighting  method. 


2.  Structure  of  TRIPIC  code 

At  first  space  meshes  are  generated  in  the  x-y  or  r-z  plane.  In 
the  space-mesh  generator  two  Poisson  type  of  equations  are  solved  [8], 
because  of  the  2.5  dimensional  code  of  TRIPIC.  The  logical  meshes  are 
composed  of  two  kinds  of  lines  with  labels  K  and  L,  shewn  in  Fig.  1 .  We 
have  to  do  the  mapping  to  the  real  space  region.  In  the  real  space 
mapped  lines  are  considered  to  be  like  equi -potential  lines.  If  the  two 
Poisson  equations  are  solved  inversely  in  the  logical  space  with  the 
appropriate  boundary  conditions,  the  numbering  of  the  space  meshes  can 
be  accomplished.  Figure  4 -a)  shews  an  example  for  the  generated  meshes. 

The  following  basic  equations  are  used  in  TRIPIC  code: 


dP/dt  =  (q/m)  (  E  +  VxB  ), 


VxE  =  -  9B/3t,  VxB  =  p (  J  +  gD/ gt  ) ,  D  =  eE 
V-D  =  p 


TRIPIC  code  employs  the  Coulomb  gauge  and,  the  Poisson  equation  for  the 
static  electric  potential  and  the  equation  of  the  magnetic  vector  one 
are  solved  on  triangular  meshes.  These  equations  can  be  solved  on  the 
cartesian  or  cylindrial  coordinate  in  TRIPIC  code.  Following  Winslow 
the  descretization  method  is  described  for  the  generalized  Poisson 
equation: 
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3<t>/3t=V(AV<i>  )  +  s 


v<fc=  t  <  <t)i-  <D  IS.*.,  -  (  4>+f  <i>  +  ■ 


The  canputational  region  covered  by  triangular  mashes  has  sets  of  six 
triangles  surrounding  a  vertex,  that  is  called  the  primary  mesh.  The 
area  of  the  secondary  mesh  is  one  third  of  the  primary  one,  shown  in 
Fig.  2.  Over  the  secondary  mesh  the  generalized  Poisson  equation  is 
integrated,  assuming  the  gradient  of  physical  quantities  constant  in 
each  triangles.  The  Gauss  law  is  adapted  to  the  first  term  at  the  right 
hand  side  of  the  above  equation.  Then  the  sum  of  the  flux  is  computed. 
Finally  we  obtain  the  descretized  equation  of  the  generalized  Poisson 
one  on  the  triangular  meshes: 


(  4>n*-  *  n)/3t  =  [  (  Jw.  <*  -  <t>  >  +  S  >/G  )9'l1/2 

”i  ”  ‘  Xi»1/2Cotei*1/5Ai-1/f=otei-l/i' 

G  =  Z  a,  £  =£  a:  area  of  the 

i=1  -6  '  secondary  mesh 

S  =  Z  s  a 

i  +  1/2  i+1/2 

By  using  this  method  basic  equations  are  also  descretized  for  the 
electromagnetic  field.  The  resulting  matrixes  of  the  basic  equations 
are  solved  by  the  direct  Gaussian  forward-  and  backward-subtraction 
method  [ 9 ] . 

By  using  electric  and  magnetic  fields  particles  are  pushed.  In 
order  to  find  the  electric  and  magnetic  fields  on  a  particle  we  should 
find  the  particle  position  in  the  logical  space.  If  the  logical  meshes 
are  numbered  in  series,  we  remember  the  position  at  the  former  time  step 
and  the  time  step  is  enough  small,  to  find  the  new  particle  position  is 
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accomplished  by  searching  the  several  meshes  just  surrounding  the  mesh 
in  which  the  particle  was  located  at  the  former  time  step.  To 
accomplish  the  interaction  between  the  particles  and  the  space  meshes 
each  particle  has  a  finite  radius  and  meshes  locating  in  this  circle 
interact  with  the  particle.  In  the  weighting  method  the  weight  is 
defined  by  the  inverse  of  the  distance  between  the  particle  and  the  mesh 
ooints(see  Fig.  6).  By  this  method  the  new  current  and  charge  densities 
are  computed. 

In  addition  TRIPIC  has  an  additional  subroutine  for  the  particle 
generation  at  the  specified  region  as  an  option.  The  generation  rule  is 
the  space-charge  limitted  condition.  For  example  electrons  and  ions 
generated  at  the  electrodes  in  an  intense-ion  diodetll ]. 

As  a  boundary  condition  the  Dirichlet-boundary  condition  is  used 
for  the  static  Poisson  equation.  For  example  at  the  anode  surface  in  an 
intense-ion  diode  the  applied  voltage  is  imposed  at  the  anode  surface. 
At  the  surface  of  the  perfect  conductor  the  condition  that  a  parallel 
electric  field  is  zero  is  imposed.  This  is  dene  even  at  the  curved 
boundary  surface  by  using  the  coupled  equation  of  the  vector  potential 
which  is  the  sunmed-up  equation  of  parallel  ccrnponent  of  equations  for 
Ax  (or  Ar)  and  Az  in  the  x-z  (or  r-z)  plane. 

As  a  summarization  of  this  section  we  present  the  time  chart  of 
TRIPIC. 


I 


n  n+1  n+2 


-x - 

X 

X 

X 

—  x- 

V 

X 

V 

X 

V 

J 

Rho 

J 

Rho 

J 

Phy 

Phy 

Es 

Es 

A 

A 

A 

B 

Et 

B 

Et 

B 

>  time  step 


Here  v  is  the  velocity,  x  the  space  coordinate,  Rho  the  charge  density, 
Phy  the  static  electric  potential,  Es  the  static  electric  field,  J  the 
current  density,  A  the  vector  potential,  B  the  magnetic  field  and  Et  the 
transeverse  electric  field. 


3.  Example  of  Computation 

In  this  section  examples  of  numerical  computations  are  presented. 

3 . 1  Mesh  generation 

The  mapping  is  done  from  the  uniform  mesh  in  the  logical  space  (K-L 
space,  see  Fig.  1)  to  the  nonuniform  real  space  mesh.  Figure  4-a)  shows 
an  example  of  mapped  meshes.  In  Fig.  1  the  computational  region  is  the 
non-hatched  region  and  the  hatched  outer  region  is  used  for  the  boundary 
treatment  of  particles  and  the  fields.  Therefore  the  mesh  size  of  most 
outer  region  has  no  need  to  be  same  with  the  mesh  size  in  the 
computational  region. 

3.2  Potential 

A  computed  electric  static  potential  in  the  region  of  Fig.  4-a)  is 
presented  in  Fig.  4-b).  The  equi -potential  lines  are  plotted  by  dotted 
lines.  In  this  case  curved  electrodes  has  constant  voltages  in  time, 
respectively  and  the  system  is  cylindrically  symmetric. 

Figure  5  shows  the  propagation  of  an  electromagnetic  single-mode 
wave.  The  propagation  of  vector  potential  is  presented  in  Fig.  5-a)  and 
Fig.  5-b)  shews  the  electric  and  magnetic  fields.  In  the  direction  of 
the  wave  vector  the  cyclic  boundary  condition  is  used. 

3.3  Particle 

In  order  to  check  TRIPIC  code  we  simulate  the  Child-Langmuir 
current  in  a  diode  gap.  Figure  7  shows  an  example  of  electron  map  in 
the  diode  gap.  In  this  case  a  gap  distance  is  0.5  cm  and  the  applied 
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voltage  is  1  volt  because  of  the  avoidance  of  the  particle  bending  by 
the  self -magnetic  field.  The  simulated  current  coincides  well  with  the 
analytical  value  which  is  obtained  only  in  the  case  of  no  magnetic 
field.  The  particle  emission  is  controled  by  the  space-charge  limit 
condition. 

4.  Conclusion 

A  particle-in-cell  code  TRIPIC  was  developed  to  do  the  particle 
simulation  in  the  irregular  boundary  by  using  the  nonuniform  triangle 
meshes.  The  electromagnetic  fields  are  solved  by  the  finite 
differencial  method  in  the  nonuniform  triangle  meshes.  The  particles 
are  pushed  by  solving  the  relativistic  equation  of  motion.  The 
interaction  between  the  particles  and  the  meshes  can  be  accomplished  by 
the  simple  weighting  method.  The  weight  is  proportional  inversely  to 
the  distance  between  the  particle  and  the  mesh.  Our  particle  has  a 
shape  of  circle  in  2  dimensional  space.  The  radius  is  comparable  with 
the  mesh  size.  The  meshes  inside  of  the  circle  can  interact  with  the 
particle.  This  weighting  method  can  be  easily  switched  to  the  other 
weighting  one. 

TRIPIC  code  can  simulate  the  space  region  with  the  irregular 
boundary.  Therefore,  for  example,  TRIPIC  is  fitted  for  the  simulation 
of  an  ion  or  electron  diode  [11]  which  has  a  curved  electrodes  to  focus 
the  beam. 

The  technique  in  TRIPIC  code  can  be  also  used  in  other  type  of  PIC 
code,  that  is  fluid  PIC  code. 
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Fig.  1  Mesh  in  Logical  Space 


Fig.  2  The  primary  and  secondary  meshes 
in  real  space 
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Fig.  5-a)  Vector  potential  in  the  case  of 
single-em-wave  propagation 
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Fig.  5-b)  Electromagnetic  wave  propagation 
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A  DOMAIN  DECOMPOSITION  AND  OVERLAPPING  TECHNIQUE 
METHOD  FOR  3D  LARGE  SCALE  NUMERICAL  SIMULATION 
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A  direct-solution  scheme  for  numerically  solving  the  3-dimensional 
Poisson’s  problem  with  arbitrarily  shaped  boundaries  V «  (AV<p)  =  S  on 
Q,  Cw+Qn-'  =  Cj  on  dQ  ,  has  been  developed  by  using  a  boundary- 
fitted  coordinate  transformation.  The  scheme  also  used  the  technique  of 
decomposing  the  closed  domain  0  into  several  hexahedron  subdomains  and  then 
overlapping  neighboring  hexahedrons  to  deal  with  complicated  geometries.  A 
large  system  of  linear  equations  deroved  from  discretizing  the  Poisson's 
equation  was  solved  by  suing  a  biconjugate  gradient  method  with  incomplete 
LU  factorization  of  the  nonsymmetric  coefficient  matrix  as 
preconditioning.  The  convergence  behavior  of  the  different  domain  decompo¬ 
sitions  was  demonstrated  for  a  numerical  experimant.  Application  to  the 
electrostatic  field  problem  in  the  electron  gun  of  a  color  picture  tube 
confirms  that  the  present  numerical  scheme  should  provide  an  efficient  and 
convenient  tool  for  solving  many  important  large-scale  engineering  pro¬ 
blems. 
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SUBTRACTION  TECHNIQUE  FOR  PLASMA  PHYSICS  AND  EVALUATING  NUM  ERICAL 

EFFECTS  IN  CODES 

Viktor  K.  Decyk 

University  of  California,  Los  Angeles 
Department  of  Physics, 

405  Hilgard  Los  Angeles,  Calif.  90024 

(  Presented  by  R.Sydora  ) 

We  have  developed  a  technique  for  simulation  of  subtle  effects  in 
plasmas.  This  technique  consists  of  performing  the  simulation  twice.  The 
first  time  one  simulates  the  background  plasma,  and  the  second  time  one 
uses  exactly  the  same  initial  conditins  plus  some  small  perturbation,  such 
as  an  extra  test  charge,  or  test  wave,  etc.  One  then  subtracts  the  results 
from  each  of  the  simulations,  to  see  only  the  effect  of  the  perturbation. 
This  method  has  been  used  to  measure  precisely  the  wake  field  excited  by  a 
single  electron  out  of  millions.  This  technique  can  be  applied  to  new 
codes  in  several  ways.  It  can  be  used  as  an  extremely  demanding  test  to 
verify  that  new  codes  are  behaving  correctly.  It  can  also  be  used  to 
explore  the  physical  behavior  of  new  models,  where  the  physical  preperties 
are  not  well  understood  beause  of  numerical  approximations  used,  such  as  in 
implicit  codes. 
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NUMERICALLY  INDUCED  STOCHASTICITY 
and 

LONG-TIME  BEHAVIOR  OF  NUMERICAL 
TRAJECTORIES  -  SMALL  AT  ANALYSIS 
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Numerically  Induced  Stochasticity 

In  a  one -dimensional  anharmonic  potential  well  tf>(x),  the  period  of  an 
orbit  is  a  function  of  its  energy.  The  true  motion  in  such  a  well  is  regular, 
since  energy  conservation  constrains  the  velocity  v  at  each  value  of  the 
coordinate  x.  Nonetheless,  when  the  orbit  is  computed  numerically,  stochastic 
behavior  can  result. 

We  consider  simple  integrators  as  mappings  from  (x,v)  at  one  time 
level  to  (x,v)  at  the  next.  With  timestep  size  A  and  <f>'  ■  8<j>/dx, 
the  leapfrog  mover  is: 

vn+1/2  -  v"'1/2  -  A  0'(xn)  ;  xn+1  -  xn  +  A  v"*1'2  . 

We  linearly  interpolate  v  to  time  level  n,  then  plot  the  pair  (xn,vn).  Equiva¬ 
lently,  we  can  write  a  variant  with  x,v  defined  at  integer  time  levels; 

v  =  v"  -  j<f>'  (xn)  ;  xn+1  =  xn  +  Av  ;  vI1+1  =  v  -  (xn+1)  . 

In  a  harmonic  well  with  natural  frequency  wQ,  the  classical  leapfrog  "dispersion 
relation"  is:  sin  wA/2  -  ±w0A/2 ,  with  stability  limit  Ac  =  2/w0  =  (period)/*-. 


We  first  noted  the  phenomenon  of  numerically  induced  stochasticity  in  the 
potential : 


For  a  >  0.0625,  this  well  has  a  central  "bump"  but  is  harmonic  with  natural 
frequency  u>0  -  42  at  large  x.  We  have  concentrated  on  the  case  a  =  1 
(fig.  1)  and  begin  with  a  description  of  the  behavior  in  this  double -well 
potential . 


The  physical  phase  plane  (fig.  2)  consists  of  closed  curves  with  one 
separatrix.  For  small  enough  A  (<  0.2),  stochastic,  behavior  is  limited  to 
orbits  near  the  separatrix  which  are  (analytically)  barely  trapped  in  one 
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The  double-well  potential  4>(x) 
obtained  when  a  =  1 . 


Fig.  2.  The  (x,v)  trajectories  of  a  set 
orbits  in  the  potential  of  fig. 


Fig.  3.  The  (x,v)  trajectories  of  a  set  of  numerically  computed 
orbits  in  the  potential  of  fig.  1;  At  =  0.25. 


sub-well,  or  are  barely  untrapped.  There  are  in  addition  large  trajectory  errors 
(without  stochasticity)  for  orbits  which  do  not  carefully  sample  the  central  bump. 


For  larger  A,  a  wide  range  of  initial  conditions  yield  a  rich  stochastic 
behavior,  with  a  complicated  island  structure.  The  situation  for  A  =  0.25  is 
depicted  in  fig.  3.  For  small  enough  energy  the  motion  is  confined  to  a  single 
sub-well  and  is  regular;  for  intermediate  values  there  is  stochasticity;  for 
large  values  there  are  large  trajectory  errors  without  stochasticity;  and  for 
very  large  energy  the  central  bump  is  negligible  and  the  motion  harmonic  (frequency 
independent  of  amplitude) .  These  effects  occur  at  timesteps  well  less  than  the 
classical  leapfrog  stability  threshold  for  the  harmonic  part  of  the  potential. 

Over  a  large  part  of  phase  space  (including  the  region  corresponding  to 
"large  trajectory  errors")  we  often  observe  a  phase-locking  phenomenon;  the 
number  of  radial  oscillations  in  the  locus  of  points  corresponding  to  an  orbit  is 
the  nearest  integer  to  the  number  of  steps  needed  to  complete  one  orbit  in  the 
harmonic  part  of  the  well.  In  the  run  depicted  in  fig.  3,  the  harmonic -well 
period  (about  4.44)  is  17.77  times  the  timestep  (0.25),  and  we  observe  18  major 
excursions.  When  A  is  0.2,  the  harmonic  period  is  22.2  steps,  and  we  observe  22 
major  excursions.  However,  while  oscillations  are  still  evident  at  A  =  0.3, 
their  number  is  44,  while  the  harmonic  period  is  14.8  steps,  or  a  third  as  long. 

Effects  of  "cantori"  are  evident  in  our  runs;  an  orbit  often  remains 
trapped  in  one  sub-region  of  a  large  chaotic  "sea"  for  a  long  time  before  it 
breaks  through  to  another  part  of  the  sea.  Runs  using  single  and  double 
precision  (Cray-1)  arithmetic  have  been  compared;  the  varying  machine  roundoff 
makes  the  number  of  steps  required  to  enter  the  various  sub-regions  of  a  large 
sea  differ,  confirming  that  the  gaps  in  the  cantori  are  indeed  "found"  by  the 
orbit  in  a  random  manner.  In  addition,  the  island  structure  is  somewhat 
different,  since  smaller  islands  can  be  resolved  in  double  precision. 

We  can  write  the  leapfrog  mover  in  a  third  form  most  suitable  for  analysis: 

xn+1  =  xn  +  Av11  -  A2/ 2  4>'  ( xn )  ;  v"+1  -  v11  -  A/2  [<*>'  (xn)  +  <6'(xn+1)] 

The  Jacobian  of  the  mapping  (xn,vn)  -♦  (xn+1,vn+1)  is  J  ■  d  (xn+1 ,  vn+1)/5  (xn ,  vn)  =  1, 
and  thus  the  leapfrog  integrator  is  an  area-preserving  map.  In  fact, 

the  leapfrog  motion  in  the  potential  well  0(x)  —  A  cos  x  is  equivalent1  to 
Chirikov's  "Standard  Map"2:  6n+1  =  +  In+1  ;  in+1  =  In  +  K  sin  9n  .  In  the 

leapfrog  scheme  as  usually  written,  we  redefine  AvI1+1/z  =  Vn+1;  then  Av11"172  =  Vn , 
and  we  obtain:  xn+1  =  xn  +  Vn+1;  Vn+1  «  Vn  +  AAzsin  xn .  Thus,  we  make  the  corres¬ 
pondence:  xo  «  ;  V  o  I  ;  AA2  o  K  ,  and  so  the  Standard  Map  theory  is 

directly  applicable  to  leapfrog  motion  in  a  sinusoidal  potential. 

The  usual  Standard  Map  orbits  (fig.  4)  are  symmetric  about  neither  0-0 
nor  1=0  (i.e.,  x=0  nor  v  =  0).  By  replacing  the  usual  leapfrog  mover  with 
our  variant  that  defines  v  at  integral  time  levels,  we  have  synthesized  a  "Sym¬ 
metrized  Standard  Map"  (fig.  5)  which  posesses  both  of  these  symmetries  and  may 
prove  useful  in  other  contexts.  In  our  experience  with  small-A  analyses  of 
various  movers  we  have  found  that  simpler  power-series  expansions  result  from 
symmetrized  schemes. 
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We  have  examined  other  wells,  and  find  that  stochasticity  is  not  limited  to 
double  wells.  For  example,  the  single-well  potential  a  —  0.0625  admits  a 
chaotic  leapfrog  orbit  for  A  —  0.75.  Chaotic  motion  is  more  pervasive  in 
double-wells;  it  sets  in  at  infinitesimal  A  because  the  phase  plane  entails  a 
separatrix  even  in  the  absence  of  finite -timestep  perturbations. 

We  have  also  examined  other  particle  moving  algorithms;  the  second  mover  we 
consider  is  an  implicit,  time-centered  mover  which  is  stable  for  all  A: 

v1*1  -  v11  -  A/2  (^'(xn+1)  +  <f>  '  ( xn ) }  ;  xn+1  -  xn  +  A/2  {v^1  +  v11}  . 

In  a  harmonic  well  with  natural  frequency  w0,  the  dispersion  relation  is: 
tan  uA/2  -  ±u0A/2.  With  a  -  1  and  A  —  0.25,  this  mover  exhibits  behavior 
similar  to  that  of  the  leapfrog.  In  the  "stiff  spring"  potential  <t>  -  x4  it 
yields  stochasticity  for  large  enough  energy  and  A. 

The  implicit  mover  is  not  (locally)  an  area-preserving  map.  The  Jacobian  of 
the  mapping  from  (x",-^)  to  (xn+1,vn+1)  is:  J  -  [l+(Az/4)<£'  '  (xn)  ]  /  [l+(A2/4)<£'  '  (xn+1)  ]  . 
Note  that  interchanging  xn  and  xn+1  amounts  to  changing  J  into  1/J,  as 
it  must  for  any  "reversible"  scheme  which  retraces  its  steps  if  the  sign  of  A 
is  reversed.  Any  change  in  phase  space  area  due  to  motion  in  one  direction  is 
made  up  for  when  the  particle  moves  in  the  other  direction.  Thus,  we  can 
expect  no  net  change  in  phase  space  area  after  integration  over  a  complete 
orbit.  For  this  scheme,  the  fact  that  J  f  1  has  no  global  consequences. 

Our  third  mover  is  a  second-order  predictor-corrector  algorithm: 

x  -  xn  +  Av”  ;  v  -  v11  -  A <y  (xn)  ; 

xn+l  _  xn  +  A(vn  +  ^)  .  v«+1  -  V**  -  |(^'(xn)  +  <6'(x))  . 

This  mover  yields  orbits  that  don't  lie  on  closed  curves  even  for  small 
timestep.  It  is  not  an  area-preserving  map,  and  does  poorly  because  the 
Jacobian  of  the  transformation  is  always  greater  than  unity:  J  ==  1  +  A 4(<£'')2/4. 

This  can't  integrate  to  zero  around  an  orbit,  and  so  the  phase-space  area  of  a 
given  set  of  initial  conditions  will  grow  continuously.  Since  only  (^'')2 
enters,  the  orbit  will  diverge  even  in  a  purely  harmonic  well,  with  a  dispersion 
relation:  w  -  w0  +  (A2w03/6)  [1  +  iAu0]  . 

The  fourth  difference  scheme3  we  consider  conserves  energy  identically; 
it  derives  from  a  symmetrized  discretization  of  Hamilton's  equations.  Let  a 
subscript  *  denote  the  advanced  time  level;  the  scheme  is,  for  any  Hamiltonian 
H(q.p). 

q*  -  q  _  1  /«Cq*.P.)  -  H(q4,p)  _  H(q,P<i)  -  H(q,p) 

A  "  2  \  P.  -  P  P.  *  P 

p»  -  p  i/H_<q*>p„)  *  H(q.p„)  ,  H(q* >p)  -  H(q,P) 
a  “  ‘  2\  q»  -  q  q.  -  q 
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Using  such  a  scheme,  the  points  of  a  numerically  computed  orbit  will  fall 
on  the  exact  phase  plane  curves  (for  a  one  dimensional  Hamiltonian) .  In  a  one¬ 
dimensional  potential  well,  the  Hamiltonian  is:  H  -  v2/2  +  $(x)  ,  and  the 
difference  scheme  becomes: 

x"*1  -  Xn  _  y"*1  +  y"  .  yn+l  -  y"  _  <ft(xn+1)  -  <j>(Xn) 

A  ”  2  ’  A  ”  xn+1  -  xn 

By  energy  conservation,  the  latter  equation  can  be  replaced  by: 

(vI1+1)2/2  -  (v")2/2  -  *(x")  -  ^(xn+1)  . 

For  most  potentials  either  form  must  be  solved  iteratively;  for  the  harmonic 
oscillator  potential  the  scheme  reduces  to  the  centered- implicit  one  described 
above.  For  <£(x)  “  x2/2,  the  orbits  are  exactly  the  circles  x2  +  v2  =  2H, 
but  there  are  phase  errors  which  grow,  in  a  single  step,  by  an  amount  =*  -A3/12. 

The  Jacobian  of  the  mapping  is:  J  -  (l-A2Q/2)/(l+AzQ„/2)  ,  where  Q  ®  (<£' -D)/(x*-x)  , 
Q„  «■  {<f>J  -D)/(x*-x)  ,  and  D  ■  (^„-^)/(x.-x) .  As  was  the  case  for  the  centered- 
implicit  scheme  above,  interchanging  x  and  x*  changes  J  into  1/J;  the  scheme 
is  reversible  and  thus  area  preserving  on  the  average. 

Other  researchers  have  advocated  the  use  of  both  energy  conserving3  and 
symplectic'*  (in  Id,  area -preserving)  schemes.  However,  at  least  in  Id  these 
desirable  attributes  seem  to  be  incompatible,  since  the  only  scheme  which 
satisfies  both  criteria  yields  a  series  of  points  along  the  true  orbit,  within 
a  uniform  (same  for  all  orbits)  rescaling  of  the  time  parameter.  Thus,  it  seems 
one  would  need  to  know  the  solution  in  order  to  devise  such  a  scheme.  The  energy 
conserving  scheme  in  Id  has  the  advantage  of  making  only  phase  errors ,  and  is 
at  least  area  preserving  on  the  average,  so  it  has  real  merit  for  at  least  a 
restricted  class  of  problems;  it  is  still  an  open  question  which  property  is 
the  more  useful  in  higher  dimensions. 

The  phenomenon  of  numerically  induced  stochasticity  has  significance  in 
several  contexts.  Firstly,  a  numerical  investigation  of  the  regions  of  phase 
space  accessible  to  an  orbit  may  lead  to  erroneous  results,  if  the  timestep  is 
too  large  or  the  mover  inappropriate.  Furthermore,  conclusions  about  orbital 
stability  based  on  numerical  integrations  may  be  erroneous,  since  neighboring 
chaotic  orbits  diverge  exponentially,  even  if  the  chaos  is  numerically  induced. 
When  studying  the  dynamics  of  a  physical  system,  one  should  demonstrate  that  any 
chaos  observed  is  not  numerically  induced.  Also,  linearized  simulations  of 
collective  phenomena  must  avoid  numerically  induced  stochasticity,  since  the 
zero-order  and  perturbed  trajectories  are  "neighboring".  Finally,  trajectory 
crossings  in  PIC  simulations  can  lead  to  enhanced  noise  and  other  errors.5 
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Long-Time  Behavior  of  Numerical  Trajectories  -  Small  At  Analysis 

One  can  construct  a  "pseudo-energy"  which  is  conserved  along  the 
numerical  trajectory.  We  define: 

*(x,v)  -  vV2  +  <£(x)  +  Alr^x.v)  +  A2<t2(x,v)  , 

and  expand  #(x.,v„)  in  a  Taylor  series  using  the  leapfrog  expressions  for 
x.  and  v,.  Defining  f  -  -<}>'  and  keeping  terms  through 

<Kx,,v„)  -  #(x,v)  +  A2[v$1'  +  f^^/Sv]  +  ... 

Thus  9  is  conserved  if  the  derivative  of  along  the  unperturbed  orbit  is  zero  : 


dir,  a*,  a* 
dF  "  '’ST  +  faF  -  °- 

This  differential  equation  has  solution  1'1  *  constant,  so  the  conserved 
quantity  is  just  the  usual  v2/2  +  <p . 

Through  the  next  order  we  find: 

»(X.,V.)  -  ’t(x.v)  +  A3[-g^  +  +  ^-]  , 

and  thus  the  equation  for  ®2  is: 
d*2/dt  -  -vff'/^  -  v3f''/12. 

With  x,v  at  interlaced  times,  the  first-order  term  A^  does  not  drop  out. 
4r1  satisfies  the  equation: 

dir1 
dF 


fMl 

2  v2f  flV 1 

lax  J 

2  l  ax2  J 

For  the  harmonic  oscillator  potential  ^(x)  -  xz/2,  we  find  that: 


tf(x,v)  -  |(  x2  +  v2  -  Axv  )  +  . . .  . 

This  amounts  to  a  skewing  of  the  constant -energy  curve  in  the  (numerical)  x,v 
plane  due  to  the  leapfrog  points  being  defined  at  "different  times  on  the  orbit". 

One  can  explicitly  solve  for  the  0(A2)  correction  to  the  conserved 
energy.  From  before,  -  -  /  dt  [  v$'<£''/^  -  v3^' ' ' /12  ].  We  use: 

dt  -  dx/v  ;  v2  -  2  (e-<£)  ;  ^ 

to  write: 

*2  -  -  J  dx  [  -  vV"/12  ] 

-  -  J  dx  [  d(*'2/8)/dx  -  €*’"/6  +  )  ■ 


-97- 


We  then  use:  d(00")/dx  -  '  +  d(0,z/2)/dx  to  write: 

n  -  ■  S  s  (  t  ■  nr  *  -  *'^2>  } 

V2^"  0 f  2 

12  24  • 

Thus,  the  conserved  quantity  is: 


f  +  4>  +  a2 


[  ] 


Along  a  numerical  orbit,  we  do  indeed  find  it  to  be  better  conserved  than  v2/2  +  0 

Using  this  "Hamiltonian,"  we  can  derive  an  equation  of  motion  which  puts 
the  particle  on  the  leapfrog  points  at  integer  times,  but  is  smooth  in  between. 
Consider  x  as  position,  v  as  canonical  momentum: 

dx/dt  —  d'S/dv  ;  dv/dt  -  -  3'f/3x . 

This  clearly  conserves  4  on  a  trajectory: 

dtf/dt  -  (30/3x)  dx/dt  +  (3V3v)  dv/dt  -  0. 

It  also  conserves  phase  space  volume;  the  volume  after  a  time  A  *»  nS  is  unchanged 
(we  take  n  infinitesimal  steps  each  of  size  6  to  get  there,  then  let  n  -+  ®) : 

x(5)  -  x  +  xS  +  .  .  .  ;  v(5)  -  v  +  <rS  +  .  .  . 


j/g-j  —  3[x(g)  ,v(3)  ]  _  1  +  53x/3x  +  .  .  .  63v/3x  +  .  .  . 

;  ~  3[x,v]  “  63x/3x  +  ...  1  +  63v/3v  +  ... 

-  1  +  5(3x/3x  +  3v/3v)  +  0(52) 

-  1  +  «<*v*  -  *„)  +  0(5Z) 

-  1  +  0(6Z) 

Thus,  after  a  time  A, 

i / * \  _  3[x(3) ,v(3) ]  3 [x(23) , v(2$) ] 

C  '  3[x,v]  X  3  [  x  ( <5 )  ,v(5)  ]  X  ••• 

-  [  1  +  0(S2)  ]  x  [  1  +  0(S2)  ]  x  ... 

There  are  n  terms,  so 

J  (A)  -  1  +  0(n52)  -  1  +  0 ( Az/n)  -  1  as  n  -►  ». 

The  equation  of  motion  along  this  "dynamical"  trajectory  is  explicitly 
expressible  (through  order  A2)  as: 


dx/dt  -  v  (1  +  Az0’ ’/6) 

dv/dt  -  -  <y  +  Az/12  (  4>'4>"  -  v*0"'  ) 
Note  that  we  have  given  up  on  dx/dt  -  v. 
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For  a  second- order  equation  for  x,  we  find: 


d2x/dt2  -  dv/dt  (1  +  aV'/6)  +  vA2/6  v<f>"  ' 

-  -  4>'  +  A2/12  (  vV ' '  -  ). 

We  compare  this  with  the  harmonic  oscillator,  for  which  the  analytic  solution  is: 

d2x/dt2  -  -  </>'  (  1  +  A2/ 12  <t>"  +  ...  )  . 

Thus,  x  and  v  will  agree,  to  order  A2,  with  the  leapfrog  solution.  The 
"dynamical"  trajectory  is  synchronized  with  the  numerical  one. 

It  is,  in  fact,  true  in  general  that  the  motion  along  the  "dynamical" 

trajectory  is  synchronized  with  the  leapfrog  motion  through  at  least  order  A3. 
Consider  the  position  as  a  function  of  time: 

x(A)  -  x  +  Ax  +  A2x/2  +  A3'x‘/6  +  ... 

-  x  +  Av(1+aV'/6)  -  aV/2  -  A3v^"/6 

-  x  +  Av  -  aV/2  +  0(A*)  . 

This  is  the  same  expression  as  that  for  leapfrog.  The  derivation  for  v  is 
similar. 
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ABSTRACT 

A  method  to  treat  particle  production  and  an¬ 
nihilation  in  paricle  code  is  proposed.  This  method 
does  not  directly  produce  particles  and  hence  can  be 
used  with  less  noise.  An  example  is  given  for  a  beat 
wave  accelerator.  The  density  change  of  0.06  %  during 
one  plasma  period  can  seriously  affect  the  excitation 
of  plasma  waves.  In  addition,  a  new  method  to 
economically  solve  atomic  processes  is  proposed,  which 
can  be  coupled  with  particle  codes.  Finally,  the  cou¬ 
pling  between  non-LTE  free  and  bound  electrons  is 
analyzed  with  a  simple  analytical  model  and  its  impor¬ 
tance  in  real  systems  is  pointed  out . 


I.  INTRODUCTION 


In  fusion  plasmas,  there  exist  many  application 
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fields  in  which  bound  electrons  inside  an  atom  is  not 
in  LTE(Local  Thermodynamic  Equilibrium).  The  non-LTE 
treatment  needs  a  solution  of  coupled  rate  equations 
and  hence  is  extremely  difficult  in  high-Z  plasmas. 
Recently,  we  have  developed  a  new  atomic  model *>  in 
which  both  charge  state  distribution  and  level  popula¬ 
tions  are  economically  calculated.  This  model  enables 
us  to  incorpolate  the  non-LTE  calculation  into  the 
hydrodynamic  processes  and  has  been  applied  to  the  X- 
ray  production  from  laser-heated  Au  plasmas2 >. 

Although  the  non-LTE  treatment  of  bound  electrons 
has  been  addressed  in  many  studies,  free  electrons 
have  been  assumed  to  be  in  LTE-Maxwell-Boltzmann 
distribution.  However,  we  find  that  there  exist  a 
wide  parameter  range  in  which  the  cross  section  of 
free-free  binary  collision  is  much  less  than  that  of 
the  electron  impact  ionization.  In  such  a  case,  even 
free  electrons  should  be  described  in  non-LTE.  There 
have  been  an  attempt  to  do  this3>,  but  the  electron 
distribution  function  was  given  outside  and  the 
modification  of  the  distribution  owing  to  the  atomic 
process  was  not  considered. 

In  this  paper,  we  discuss  the  possibility  to  es¬ 
tablish  a  new  particle  model  to  describe  non-LTE  free 
electrons  incorpolating  interactions  with  non-LTE 
bound  electrons  :  ionization  and  recombination.  In 
section  II,  we  propose  a  simple  algorithm  to  treat  an- 
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nihilation  and  production  of  particles. 


In  section 


III,  this  scheme  coupled  with  a  relativistic 
electromagnetic  particle  code,  which  is  essentially 
the  same  as  ZOHAR4>,  is  applied  to  a  beat  wave 
accelerator5  >  6  > .  In  section  IV,  the  review  of  recent 
atomic  models  are  given.  In  section  V,  the  coupling 
between  the  atomic  process  and  the  particle  distribu¬ 
tion  is  treated  by  a  simple  model  and  the  model  is  ap¬ 
plied  to  a  problem  on  heat  flux  in  laser-produced 
plasmas . 


II.  PARTICLE  CODES  WITH  ANNIHILATION  AND  PRODUCTION 

Particle  codes  have  been  used  in  many  areas  ; 
electromagnetic  particles  in  plasma  physics,  gas  par¬ 
ticle  in  astrophysics,  fluid  particles  in  fluid 
mechanics  and  so  forth.  The  code  can  be  easily  imple¬ 
mented  and  is  tolerable  in  large  distortion  or 
disturbance.  In  these  applications,  however,  there 
have  been  little  attempt  to  incorpolate  the  annihila¬ 
tion  and  production  of  particles  in  the  codes.  This 
is  because  the  annihilation  and  production  of  par¬ 
ticles  require  inefficient  and  tedious  caculations  and 
sometimes  cause  a  serious  noise.  In  this  section,  we 
propose  a  simple  algorithm  to  incorpolate  the  above 
process  into  particle  codes. 

In  most  cases,  a  particle  is  produced  at  some 
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point  having  a  velocity  distribution.  Let  us  define 
the  production  rate  of  this  disctribution  by 
df(x,v)/dt.  In  order  to  replicate  this  production,  a 
number  of  particles  must  be  generated  ;  velocity  dis¬ 
tribution  must  be  replicated.  Furthermore,  if  df/dt 
is  quite  small,  a  particle  will  be  generated  very 
rarely  or  discontinuosly  leading  to  noise  because  mass 
of  generated  particle  cannot  be  small  ;  if  the  par¬ 
ticle  mass  is  taken  to  be  small,  number  of  particles 
generated  will  be  extremely  large.  In  the  present 
scheme,  the  particle  is  not  directly  generated.  In¬ 
stead  of  generating  directly  the  paricles,  the 
produced  particles  are  stored  into  a  cell  as  in  a  form 
of  a  velocity  distribution  f i (v) ,  which  is  advanced  in 
time  as 

fiDtl(v)=fi"(v)+df(xi ,v)/dt.A  t,  ( I I — 1 ) 

where  x,  is  the  location  of  the  i-th  cell  center  and  n 
du..otes  the  time  step,  and  A  t  is  the  increment  of  the 
time  step.  The  distribution  does  not  propagate  in 
space  and  has  no  response  to  forces  generated  or 
imposed.  In  our  algorithm,  new  particles  will  not  be 
produced  unless  the  severe  condition  described  later 
is  met.  Instead  of  generating  new  particles,  the  mass 
of  existing  particles  which  are  in  the  i-th  cell  are 
increased  by  the  following  procedure; 

Np“*“(xp,vp)rNp® ii(xp,vp)+a  f iB,1(vp)A  v,  ( I 1-2 ) 
where  Xp  and  vP  are  the  location  and  velocity  of  the 
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particle,  and  A  v  is  the  velocity  volume  ;  fi  is 
divided  into  many  cells  in  velocity  space.  In  turn, 
the  stored  particle  decreases  according  to 
f  »■  — (vp)  =  (l-a  )  f  i  “  *  1  ( vP  )  ,  (H-3) 

in  one  time  step.  In  the  example  given  in  section  III, 
a  is  set  to  1  or  0  (see  also  Eq.(II-4)). 

We  should  impose  the  maximum  and  minimum  size, 
Nali  and  Naia  on  the  particles.  Naax  avoids  too  large 
particle  which  may  become  a  source  of  noise.  Then 
a  =0  if  ( 1 1  —  4  ) 

Once  the  mass  increase  stops  as  in  Eq.(II-4),  fi(vp) 
corresponding  to  this  particle  may  increase 
intolerably.  This  situation  is  unavoidable,  par¬ 
ticularly  for  cold  streaming  electrons  in  uniform  ion 
background  ;  some  of  these  beam  electrons  collide  with 
bound  electrons  inside  ions  and  a  new  "non-drifting” 
electron  should  appear.  Since  there  is  no  electron  in 
this  velocity  range,  a  new  particle  must  be  generated. 
However,  if  f i  is  used,  this  procedure  can  be  simply 
performed.  Thus,  if  fi(vp)A  Nm,  a  particle 
having  vp  and  Maia  is  generated  inside  the  i-th  cell. 

Annihlation  of  particles  is  very  simple  in  this 
procedure.  It  should  be  noticed  that  a  particle  will 
not  disappear  even  if  Np<Nn  D , 

The  only  problem  is  how  to  store  this  distribu¬ 
tion  with  less  computation  memory  ;  space  times 
velocity  dimensions  is  required. 
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III.  APPLICATION  TO  BEAT  WAVE  ACCELERATOR 


In  this  section,  we  give  a  very  simple  applica¬ 
tion  of  the  procedure  given  in  the  previous  section. 
That  is  a  beat  wave  accelerator5  > 1 > .  There  have  been 
a  number  of  particle  simulations  on  this  problem, 
while  none  has  incorpolated  ionization  effects  into 
the  particle  code.  In  some  of  experiments  on  the  beat 
wave  accelerator,  medium-Z  materials  are  used  as  a 
plasma  source.  Even  in  this  medium-Z  plasmas,  ions 
may  not  be  fully  ionized.  During  the  heating  of 
plasmas,  the  ion  density  gradually  increases  in  time. 
Although  the  rate  of  density  change  is  very  small,  it 
can  affect  the  excitation  of  plasma  waves.  This  was 
first  pointed  out  by  J.P. Matte  et  al . 7 > . 

It  is  well  known  that  plasma  waves  excited  by 
beat  wave  saturate  because  of  a  relativistic 
detuning6’.  This  is  due  to  the  increase  of  y  =[1- 
(v/c)2]-,/2,  where  v  and  c  are  the  speeds  of  a  par¬ 
ticle  and  light,  respectively.  The  idea  by  J.P. Matte 
et  al .  is  to  compensate  this  reduction  of  plasma 
frequency  by  increasing  the  ion  density  through 
ionization.  This  change  in  $  u  p  is  approximately 
written®  >  as 

6  u>  */(*>  p  =  -l/2  (9a,*; '8)  2'3  (III-l) 

where  ai ,  aj  (<<1)  are  the  normalized  quiver 
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velocities  eEi/(mu  ic),  eEi/(mu  2c).  Hence  the  den¬ 
sity  change  <5  n  which  compensates  for  the  change  given 
by  Eq . ( I I I - 1 )  should  be 

6  n/n  =  (9aia2/8)*/s  (I1I-2) 

Since  the  time  to  saturate  is 

co  Pr  s  =(4/31/3)  (aiaj/4)  ~2/3>  <  I X I  —  3  ) 

the  rate  of  density  change  is 

i  n/nw  pt  ,:(3'/3/4)(9/32)  *  /  3  ( a  1  a  2  )  *  /  a  .  (HI-4) 

For  CO2  laser  (10.6  and  9.6  n  m)  at  1014W/cm2,  aia:  = 
0.0074  and  hence  Eq.(III-2)  gives  8  n/n  =  0.04.  This 
indicates  that  this  sma1 1  change  of  density  can  affect 
the  excitation  of  plasma  waves. 

Since  o>  Pr  »  =  184  and  8  n/nta  Pr  s  =  2.2xl0-4,  the 
density  changes  by  only  0.022%  in  one  plasma  period  ( co 
P *  1  ) .  However,  even  this  small  change  can  be  success¬ 
fully  described  with  the  algorithm  given  in  section  II 
Let  us  give  an  example  ;  ai=a2=0.1,  co  i/w  p=10.9, 
<0  2/co  P  =  9.9,  the  system  size  is  6.35c/cu  P,  total  spa¬ 
tial  meshes  256.  The  initial  electron  temperature  is 
m*c2/25.  The  ejected  electrons  through  ionization  are 
stored  into  fi(v)  as  in  Eq.(II-l)  ;  the  velocity  space 
is  divided  inot  10x10  meshes  for  the  range  -1^  vx/c, 
Vj/ci  1.  N ■  a  1  and  Nn  B  are  set  to  2.0No  and  0.5No,  No 
being  the  initial  particle  mass. 

Figure  1 1 1 —  1  shows  the  time  evolution  of  par¬ 
ticles  in  phase  space  (x,ux)  ;  this  is  the  case 


without  ionization. 


Figure  III-2  shows  the  time 


evolution  of  the  maximum  electrostatic  potential  in¬ 
duced  by  beat  wave  ;  III-2a  and  2b  correspond  to  the 
cases  without  and  with  ionization,  respectively.  In 
the  latter  case,  ionization  rate  n/noi  Pr  is  6.25  x 
10'4,  whereas  the  value  given  by  Eq.(III-4)  is  3.33  x 
10-4.  As  clearly  seen  from  these  figures,  even  this 
small  change  of  the  density  can  seriously  modify  the 
beat  wave  acceleration .  This  may  also  be  true  for  the 
induced  Raman  process  in  laser-produced  plasmas. 


IV.  ATOMIC  PROCESS 


In  this  section,  the  method  to  solve  atomic 
processes  is  briefly  reviewed.  The  basic  equations 
for  those  processes  are  the  rate  equation  for  the 
population  Ni ,B  of  an  ion  in  charge  state  Z  with  an 
excited  electron  in  the  n-th  level; 


"  *  »'  *' 
1  >  i  / 


This  equation  is  not  economical  when  it  is  coupled  to 
other  equations  which  determine  the  plasma  parameters. 

In  practice,  it  is  better  to  use  an  approxima¬ 
tion  to  this  equation  and  to  make  it  realistic.  There 
exist  three  approaches  at  hand.  The  first  one  is  the 
average  ion  model9' 10 >,  where  Nz,n  is  averaged  over  Z 
and  then  it  is  sufficient  to  solve  only  the  equations 
for  the  level  populations  of  an  average  atom.  The 
equation  to  be  solved  in  this  case  is 


diNP.)/dt=R'N.ffQ,-I.f/.ffP,+ffl  2  AmPmQ,~  2  Am,P,Qm  j  +  2  C^N.NPmQ, 

l  m  >n  m' <*  \  m  >n 

-  2  cS,;N.NP,Qm  +  2  C^N,NPmQ,~  2  C^.N,NPnQm. ,  (VU~  2  ) 

m’<*  m  <■  m'>n  ‘  ' 


The  derivation  of  this  equation  is  described  in  Ref.l 
in  detail . 

As  an  alternative  approach,  if  a  suitable  assump¬ 
tion  for  the  populations  of  the  excited  level  can  be 


made,  the  equations  required  to  be  solved  are  those 
for  N*  only;  N*  being  the  abundance  of  the  ion  Z.  The 
equation  for  this  is 

dN,/dt=-IINtN,-RlNtN,+Jx_iN,_iNt 

( IV-3 ) 

+R,  +  iN,+iN.  > 

In  this  direction,  Salzmann  and  Krumbein11*  proposed  a 
form  of  the  excited-level  population  such  as 

N x  ,  i/N i  =  NoAexp(-E»  ,  B/kT)  for  excited  levels 
N * , »/N i  =  No  for  a  ground  level 

where  A  is  so  chosen  that  the  population  becomes  close 
to  the  value  calculated  by  a  detailed  rate  equation. 
Here,  No  is  determined  from  2  »Ni,«/Nt=l  and  E»lB  is 
the  level  energy.  However,  this  is  not  always  pos¬ 
sible  for  high-Z  plasmas,  since  we  have  no  method  to 
determine  the  adjustable  parameter  A.  Busquet 

proposed  a  mixed  model  where  only  relatively  lower  ex¬ 
cited  levels  are  calculated  by  the  rate  equation 
whereas  highly  excited  levels  are  assumed  to  be  of  the 
Boltzmann  type.  However,  he  needed  a  further  ap¬ 
proximation  in  the  coronal  limit. 

The  third  approach  is  the  hybrid  atom  model *> 
developed  recently,  which  is  a  combination  of  the 
former  two  methods.  The  basic  principle  of  this  ap¬ 
proximation  comes  from  the  fact  that  transitions  be¬ 
tween  bound  states  proceed  faster  than  those  between 
free  and  bound  states.  Hence,  the  level  populations 


of  the  excited  levels  can  be  approximately  estimated 
by  quasisteady  equations  of  bound-bound  transitions 
for  an  ion  of  charge  Z: 

+C£JV.)/(0V,>  ( IV_4  ) 
=(gm/gm)F(Em,Em,N.)  , 

Equation  (IV-4)  shows  that  the  normalized  population 
Ni.a/ga  can  be  estimated  if  only  the  level  energy  Ea 
is  known.  Although  there  exist  transition  processes 
which  may  destroy  the  above  characteristics,  Eq.(IV-4) 
motivates  us  to  use  an  energy-functional  form  of  the 
normalized  population  Nt,a/gt,n  :  if  N2,a/g2,a  is  only 
a  function  of  level  energies  and  does  not  depend  ex¬ 
plicitly  on  Z,  the  normalized  populations  constructed 
from  the  average  ion  can  be  rescaled  and  used  for 
other  ions  in  a  different  charge  state.  It  should  be 
noticed  that  we  do  not  use  Eq.(IV-4)  but  use  Eq.(IV-2) 
in  the  hybrid  atome  model. 

In  order  to  clarify  the  idea  in  detail,  let  us 
examine  the  level  dynamics  of  various  ions  by  solving 
the  full  rate  equation  (IV-1)  for  an  aluminum  plasma  ; 
it  is  possible  to  solve  (IV-1)  in  such  a  low-Z  plasma. 
In  Fig. IV-1,  the  reduced  population  probability 
Wj.i/gj,.  (=  2  «WZ|B=1)  is  plotted  versus 
the  excitation  energy,  which  means  the  level  energy 
measured  from  the  ground-state  level.  For  Ni=10*° 
cm-3,  the  probability  is  close  to  the  Boltzmann  type 
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Wz , B=gz , Bexp(-E„/kT* ) .  Although  for  Ni=10»«  cm-3  it 
is  also  close  to  the  Boltzmann  type,  the  temperature 
of  the  bound  electrons  is  much  lower  than  the  free 
electron  temperature  T«.  In  each  case,  the  reduced 
probability  for  various  ions  of  +3  to  +10  can  be  well 
characterized  universally  by  one  curve,  although  this 
curve  depends  on  the  physical  parameters  as  shown  by 
two  curves  in  Fig.(IV-la). 

If  this  characteristic  is  taken  into  account,  the 
common  behavior  appearing  in  Fig.(IV-la)  can  be  satis- 
>  factorily  described  by  the  average-ion  model. 

Actually,  the  reduced  electron  population  probability 
defined  later  shows  the  same  characteristics  as  shown 
in  Fig.(IV-lb).  Accordingly,  once  the  electron 

population  P„  of  the  average  ion  is  obtained  from 
Eq . ( IV-2 ) ,  the  function  Y»(E)  (the  reduced  electron- 
population  probability)  is  constructed  from  P„  as  fol¬ 
lows  ;  at  some  discrete  points 

Y  a { E  n  * )  =  P»/gi2  nPn  ( n= 1 , 2 , — n . a  ,  )  (IV-5) 

and  in  other  regions  Ya(E)  is  exponentially  interpo¬ 
lated  from  Eq.(IV-5)  :  this  exponential  interpolation 

is  justified  by  the  characteristics  appearing  in 
Fig.(IV-l).  In  Eq.(IV-5)  EB*  is  the  excitation  energy 
of  the  average  ion  and  the  summation  is  taken  over  the 
ionizing  shells.  This  energy-functional  form  can  be 
used  to  generate  the  population  probability  even  for 
the  ion  in  a  charge  state  Z  different  from  Z»  of  the 


k 
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average  ion  : 


W’ 2  ,  i-Ya(E)gz i ■ .  ( IV-6  ) 

It  should  be  noticed  that  the  physical  meaning  of 
W'z,n  in  Eq.(IV-6)  differs  from  Wz,n  in  Fig.(IV-la)  : 
the  former  is  derived  from  the  electron  population  but 
the  latter  from  the  ion  population  Nz,n.  However, 
both  characteristics  are  similar  as  seen  from 
Figs.(IV-la)  and  (IV-lb).  This  can  be  justified  by 
remembering  the  nature  of  the  average  ion  model  and  CR 
( Coll isional-Radiative )  model.  In  the  CR  model,  only 
one  electron  is  excited  in  some  level.  Accordingly, 
if  ions  in  various  charge  states  are  collected 
together  into  a  fictitious  ion,  the  electron  popula¬ 
tion  PB  in  the  n-th  level  of  this  ion  will  reflect  the 
ion  population  of  the  original  "real"  ions  Nz,n  with 
an  electron  in  the  n-th  level. 

In  the  next  step,  the  ionization  and  recombina¬ 
tion  rate  coefficients  are  summed  up  with  respect  to  n 
using  N z  and  W!in:W’!in,  Then  Eq.(IV-3)  for  Nz  is 
solved.  With  this  model,  all  regions  including  LTE , 
corona,  and  those  intermediate  to  them  can  be 
described  without  any  adjustable  parameter. 

Since  Eq.(IV-3)  is  in  a  tridiagonal  form,  the 
numerical  procedure  to  solve  it  is  quite  easy.  In  the 
present  model  the  main  computation  time  becomes  the 
time  to  solve  only  Eq.(IV-2)  and  is  computationally  as 
fast  as  the  average-ion  model.  If  the  ionic  charge 


1 


becomes  large,  then  the  validity  of  the  present  model 
is  much  improved  because  the  range  of  Z  that  satisfies 
E  =!Z-Z»!/Z»  <<  1  becomes  wider  and  the  fractional 
energy  change  (A  E/E)  for  these  ions  is  on  the  order  of 
E  ;  this  small  change  of  energy  justifies  the  inter¬ 
polation  in  constructing  the  energy  funcion  Y». 

In  order  to  justify  the  model,  let  us  give  an 
example.  Figure  (IV-2)  shows  the  relative  abundance 
of  the  charge  state  of  aluminum  at  a  steady  state  for 
Ni  =  10 20  cm-3  and  various  electron  temperatures. 
Here,  the  rate  equation  [Eq.(IV-3)]  for  N*  includes 
radiative,  three-body,  and  dielectric  recombinations, 
and  collsional  ionization.  In  addition,  collisional 
excitations  and  deexcitations ,  and  radiative 
deexcitations,  are  also  included  in  the  average-ion 
model  [Eq. ( IV-2 ) ] .  The  results  predicted  by  the 
hybrid  atom  model  [Eqs.(IV-2)  and  ( IV-3 ) ] ( sol id  line) 
agree  quite  well  with  those  given  by  Duston  et  al . 1 3 > 
(CR  model  shown  by  the  dashed  line)  for  ZJi  10.  In 
Ref. 13,  only  levels  for  Z£  10  were  calculated  by  the 
rate  equation,  whereas  for  Z<10  ions  only  the  ground- 
state  levels  were  taken  into  account.  The  small  dif¬ 
ference  (about  40%  at  most)  between  the  CR  model  and 
our  hybrid-atom  model  is  within  difference  caused  by 
the  different  rate  coefficients. 

The  superiority  of  the  hybrid  atom  over  the 
average  atom  appears  in  the  calculation  of  radiation 
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transport  in  high-Z  plasmas.  Figure  IV-3  illustrates 
the  main  feature  of  this  process.  Let  us  consider  x- 
ray  emission  at  some  spatial  point  xi  and  absorption 
at  x 2 •  In  the  average-ion  model,  only  one  fictitious 
ion  having  Z*  exists  at  each  space  point,  so  that  the 
level  energy  and  hence  the  x-ray  spectral  energy 
changes  from  space  point  to  space  point  depending  on 
the  average  charge  Z*  which  follows  physical 
parameters,  the  density  and  the  temperature. 
Consequently,  the  line  radiation  emitted  in  one  region 
can  not  be  absorbed  by  the  same  spectral  line  in 
another  region  as  shown  in  Fig. IV-3 (a).  In  actual 
plasmas,  however,  there  exists  an  ion  having  the  same 
charge  Z  in  all  regions  although  the  abundance  Nz  may 
change  from  space  point  to  space  point.  This  means 
that  a  spectral  line  emitted  in  one  region  is  always 
absorbed  by  the  same  line  in  any  other  region  of  space 
as  shown  in  Fig. IV-3 (b).  The  hybrid-atom  model  can 
describe  this  process. 

The  line  shift  due  to  change  of  Z  is  significant 
in  high-Z  material  because  the  energy  difference  be- 


tween 

Z  and 

Z+l 

ions  is 

about 

21  B 

Z/n2 

and 

hence 

is 

about 

50  eV 

for 

Z=30  and 

II 

c 

I. 

being 

13.6 

eV. 

This 

value 

is  far 

beyond  the 

line 

width 

broadening . 

The 

result‘>  with  the  hybrid-atom  model  indicates  that  the 
charge  state  of  Au  plasmas  distributes  over  a  rela¬ 
tively  large  width  in  Z  ( A  Z=5).  This  result  also 
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agrees  with  the  experimental  data14*.  The  real  lines 
may  spread  over  300-400  eV  (n=4)  in  spectral  range  ; 
the  line  of  the  average  ion  may  locate  at  the  center 
of  this  group  of  lines. 

As  a  simple  example,  we  solve  line  transfer  equa¬ 
tion  as  well  as  the  atomic  process.  We  use  exponen¬ 
tial  profiles  for  the  density  and  temperature  in  space 
making  the  pressure  uniform  :  the  density  and  tempera¬ 
ture  ranges  from  16x10*°  to  1x10*°  cm**  and  100  to 
1600  eV.  The  scale  length  of  these  plasmas  is  10  n  m. 
Figures  IV-4(a)  and  4(b)  are  obtained  from  the 
average-ion  model  and  the  hybrid-atom  model, 
respectively.  In  addition  to  the  atomic  processes 
given  above,  photoexcitations  and  radiation  transport 
are  included1).  In  the  figure,  the  solid  curves  are 
the  results  in  which  the  energy  space  is  divided  into 
groups  of  5  eV  width  and  line  broadenings  by  various 
processes  are  neglected  for  simplicity  :  one  curve  is 
the  direct  data  and  another  is  further  averaged  over 
50  eV  for  an  illustrative  purpose.  In  Fig.IV-4(a),  in 
addition  to  the  above  curves,  we  draw  the  dashed  curve 
that  is  the  result  using  the  groups  of  50  eV  width 
broadened  artificially.  The  comparison  between  the 
solid  curves  in  Fig.IV-4(a)  and  4(b)  shows  that  the 
emission  intensities  and  hence  opacities  with  the 
hybrid-atom  model  are  larger  than  those  with  the 
average-ion  model.  The  difference  is  concluded  to 
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arise  from  the  charge-state  distribution  illustrated 
in  Fig.IV-3.  It  is  possible  to  approximately  incorpo- 
late  this  effect  into  the  average-ion  model.  If  the 
line  width  is  artificially  broadened  over  the  range 
determined  from  the  charge  state  distribution  as  in 
Fig.IV-4(a)  (dashed  curve),  the  result  comes  closer  to 
that  with  the  hybrid-atom  model. 
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V.  COUPLING  BETWEEN  NON-LTE  FREE  AND  BOUND  ELECTRONS 


In  this  section,  we  will  describe  an  example  of 
coupling  between  free  and  bound  electrons  both  in  non- 
LTE  state.  In  order  to  clarify  the  physics,  the  model 
is  largely  simplified.  However,  its  basic  principle 
can  be  used  in  the  particle  code  described  in  section 
II . 

In  the  formulation  of  section  IV  and  Ref.l,  free 
electrons  are  assumed  to  be  in  LTE.  This  assumption 
appears  in  the  rate  coefficient  R  as  follows; 

R=<  <r  v>=  vfv*dv,  (V-l) 

where  <r  is  the  cross-section  and  f  the  velocity  dis¬ 
tribution  function  of  free  electrons  ;  f  is  normalized 
to  satisfy 

^fv*dv  =  1,  ( V-2 ) 

and  is  taken  to  be  the  Maxwell  distribution  in  LTE. 
In  this  section,  however,  we  determine  this  distribu¬ 
tion  function. 

It  is  easy  to  derive  the  rate  of  change  in  the 
number  of  free  electrons  owing  to  radiative  capture  ; 
the  number  of  free  electrons  captured  within  a 
velocity  width  A  v  at  v  is  given  by  ( <5  R/<?  v)A  vNNe  = 
<T  *fv3A  vNNe.  On  the  contrary,  the  expression  for  the 
impact  ionization  is  more  complicated  because  there 
are  three  terms  in  this  process  ;  1)  velocity  change 


\ 


from  v  to  v*  of  a  projectile  electron  ,  2)  from  v'  to 

v  of  the  same  electron,  and  3)  the  birth  of  a  free 
electron  with  v  (ionization).  If  we  denote  the  ini¬ 
tial  and  final  velocities  of  the  free  electron,  and  an 
emerging  electron  velocity  by  vj,  vt,  and  vt, 
respectively,  the  energy  conservation  requires 
vi*  =  2I/m.  +  Vf*  +  vb*,  (V-3) 

where  I  is  the  ionization  energy.  Since  the  cross 
section  of  binary  collision  between  free  electrons  is 
much  larger  than  that  of  impact  ionization  at  low 
energy  (close  to  ionization  energy  I,  see  Fig.V-1), 
we  need  not  to  take  into  account  the  third  process 
(birth  of  an  electron)  ;  the  ejected  electron  may 
quickly  thermalize  through  binary  collision  which  can 
be  described  as  -v  (f-fs),  where  f«  is  the  equilibrium 
Maxwell  distribution.  y  is  given  by  a  «evN*  where  a 
e «  is  the  cross  section  of  free-free  collision. 

Furthermore,  we  can  also  assume  I/me  <<  Vi!,  Vf*. 
Thus,  we  get  an  approximate  equation  ; 

df/dt  =  -ff  RvNef(v)+(T  i  v  '  N  «  f  ( v '  )  -  <y  ivN*f(v) 

-V  (v) (f-fn)  ( V-4 ) 

where  v*  is  the  initial  velocity  of  a  projectile 
election.  If  we  set  v’=v  +  <A  v>  and  expand  Eq.(V-4) 
in  <A  v>  <<  v,  v’ ,  then  we  get 

df/dt  =  -<y  nvNef+tx  i<Av>N*  9(fv)/<5v 

-<I  •  »vN *  (  f - f  m  )  .  ( V-5  ) 

where  <A  v>  can  be  approximately  given  by  I/m*v. 
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As  inferred  from  the  relative  magnitude  of  <y  in 
Fig.V-1,  the  electrons  having  high  energy  may  be  much 
affected  by  ionization  a  i.  Figure  V-2  shows  the 
nearly  quasi-steady  solution  of  Eq.(V-5)  for  Au  plas¬ 
mas  at  T=lkeV,  Z*=41.  This  change  of  velocity  dis¬ 
tribution  is  quite  serious  in  calculating  the  thermal 
flux  and  other  transport  coefficients. 
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FIGURE  CAPTIONS 

Fig.III-1  :  Phase  space  plot  of  particles  at  u  pt  = 
O(top),  lOO(middle),  and  200(bottom). 

Fig.III-2  :  Time  evolution  of  electrostatic  potential 
(a)  without  and  (b)  with  ionization. 

Fig .  IV- 1  .'  The  reduced  population  probability  vs  the 

excitation  energy  for  T=50  eV  and  (i)  N=1020  cm*3  and 
(ii)  N=1019  cm*3.  (a)  The  reduced  ion  population 

probability  vs  the  excitation  energy  with  CR  model  for 
different  charge  states  of  ions.  (b)  The  reduced 

electron  population  probability  vs  the  excitation 
energy  with  the  average-ion  model. 

Fig.IV-2  :  The  dependence  of  aluminum  abundance  in 
various  charge  states  on  electron  temperature  for  an 
ion  density  of  1020  cm*3.  The  solid  and  dashed  lines 
are  the  results  from  the  hybrid-atom  model  and 
collisional-radiative  model,  respectively. 

Fig.IV-3  :  The  difference  between  the  radiation 
transport  in  (a)  average-ion  model  and  (b)  hybrid-atom 


model . 


Fig.IV-4  :  The  emitted  x-ray  intensity  vs  photon 
energy  from  a  typical  laser-produced  Au  plasma. 

(a)  average-ion  model,  (b)  hybrid-atom  model. 

Fig.V-1  :  Comparison  between  free-free  and  bound- 
free  collisions  ;  here  T=I  is  assumed.  When  T<I, 
free-free  cross  section  becomes  steeper  and  shifts 
towards  left.  £  is  the  number  of  bound  electrons. 

Fig.V-2  :  Quasi-steady  solution  of  Eq.(V-5)  for  Au 
plasmas  at  T=1  keV,  Z*=41.  Dashed  line  shows  the  Max¬ 
well  distribution  at  T=1  keV. 
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A  long-standing  goal  in  plasma  simulation  has  been  a  method  which  could 
treat  both  detailed  kinetic  physics  and  smooth  large-scale  physics  in  an 
efficient  and  natural  way.  Until  recently,  particle  simulations  were  applied 
almost  exclusively  to  problems  of  "microscopic"  physics,  where  only  a  small 
part  of  the  plasma  was  modeled.  With  the  advent  of  implicit  particle 
simulation  techniques,1'5  one  can  now  treat  systems  of  many  Debye  lengths  and 
follow  them  for  many  plasma  oscillation  times.  However,  the  price  paid  for 
this  capability  is  a  restriction  on  the  allowed  spatial  resolution  (an  accuracy 
constraint).  Thus,  at  the  current  state-of-the-art  (e.g.,  in  codes  such  as 
TESS6  or  AVANTI7)  a  small  timestep  is  still  necessary  whenever  the  system 
incorporates  a  physically- important  small  spatial  or  temporal  scale  anywhere 
within  its  domain. 

We  are  working  to  develop  a  new  particle- in-cell  plasma  simulation 
technique  which  relaxes  these  restrictions  and  would  be  suitable  for  strongly 
inhomogeneous  problems  involving  a  wide  range  of  space  and  time  scales.  The 
plasma  in  any  part  of  phase  space  would  be  advanced  on  its  own  natural  scales . 
Of  course,  it  seems  only  natural  to  advance  each  particle  with  its  own, 
independent,  series  of  timesteps;  one  could  imagine  using  (for  example)  an  ODE 
solver  such  as  LSODE.  The  major  difficulty  in  doing  this  has  been  the 
necessity  of  processing  the  particles  in  synchrony  due  to  the  requirement  of  a 
self-consistent  field.  We  have  developed  in  outline  an  algorithm  which  may 
overcome  this  difficulty.  The  code  would  advance  the  particles  in  blocks  k, 
each  with  an  associated  timestep  At^.  For  each  region  of  phase  space,  the 

nominal  timestep  (and  possibly  the  mesh  spacing)  could  be  chosen  in  an  adaptive 
manner.  A  major  improvement  in  economy  comes  about  because  the  majority  of  the 
particles  are  not  processed  during  any  given  step;  for  suitable  problems  this 
gain  may  be  two  orders  of  magnitude. 

There  are  many  areas  of  plasma  physics  where  such  a  capability  would  be 
highly  desirable ,  including  sheaths  at  the  interface  of  a  plasma  with  a  probe 
or  wall,  collisionless  shocks,  double  layers,  and  a  wide  variety  of  astro - 
physical  problems.  Even  in  a  relatively  tractable  problem  such  as  the  bump-on- 
tail  instability,  significant  gains  might  be  realized  by  pushing  only  the  fast 
particles  with  a  small  timestep;  the  more  numerous  bulk  particles  would  be 
advanced  less  frequently.  Similar  methods  might  prove  applicable  to  problems 
of  gravitating  systems  and  to  particle- in-cell  fluid  modeling. 

In  the  scheme  we  are  testing,  particles  are  advanced  only  every  j  =  2m 
steps  (for  some  integer  m  >  0)  If  they  reside  in  a  region  of  phase  space  where 
there  are  no  large-amplitude  short -wavelength  rields  (kvAt  and  wtrapAt  are 
moderate) .  In  such  a  region  we  may  additionally  employ  a  coarse  mesh  or  strong 
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spatial  filtering.  The  large  timestep  used  in  such  a  region  should  help  avoid 
the  finite-grid  instability  common  to  explicit  coder  with  coarse  meshes. 


The  direct- implicit  field  equation  is:  V*(l+x)V$  “  P.  with  x(x)  =  fw^At2 
and  p(x)  a  free -streaming  density.  It  is  solved  over  the  entire  domain  every 
step.  In  this  way,  the  deposition  of  charge  occurs  implicitly ,  one  step  earlier 
than  in  an  explicit  code.  We  allow  a  block  to  deposit  its  information  j  time 
levels  ahead  of  the  current  one;  this  information  is  then  interpolated  backward 
in  time  to  yield  the  data  needed  to  produce  the  field  a  single  time  level  ahead. 

As  particles  move  about  in  phase  space,  it  is  necessary  to  change  their 
timesteps  (move  them  from  block  to  block) .  We  allow  changes  in  At  by  no  more 
than  a  factor  of  2,  no  matter  where  the  particle  is;  we'll  catch  up  on  the  next 
step  anyhow,  and  the  logic  is  simpler  this  way.  To  facilitate  doubling  or 
halving  the  step  size  "between  steps",  we  employ  a  variant  of  the  "dl"  implicit 
particle  advance  with  all  key  quantities  defined  at  integral  (not  staggered) 
time  levels.  This  should  allow  us  to  preserve  second  order  accuracy  in  time. 

A  typical  grid  might  look  like: 


(for  electrons)  At  -  45t  (nominally)  At  -  25t  At  =  St 

j  -  4  j  -  2  j  =  1 

For  each  value  of  j  there  are  j  blocks: 


Block  el:  Electrons:  push  every  step 

Block  e2:  Electrons:  push  on  even-numbered  steps 

Block  e3:  Electrons:  push  on  odd-nvimbered  steps 

Block  e4:  Electrons:  push  on  steps  with  (step  number  mod  4)  -  0 

Block  e5:  Electrons:  push  on  steps  with  (step  number  mod  4)  -  1 

Block  e6:  Electrons:  push  on  steps  with  (step  number  mod  4)  =  2 

Block  el:  Electrons:  push  on  steps  with  (step  number  mod  4)  -  3 


Blocks  il-i7  are  similar,  but  contain  ions. 

Let  us  abbreviate  "time-level"  by  "tl;"  then  the  code  at  timestep  7 
processes  the  active  blocks  (dots  denote  interpolation  in  time  of  p  and  x)  '■ 

Time  level  --3  4  5  6  7  8  9  10  11 


Blocks  1  /x~\  /p^x\ 

3  /  w  \  /p,x  f...\ 

7  /  x7^  \  /p,x  7 . \ 
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The  other  blocks  were  advanced  on  earlier  steps ,  and  we  need  only  interpolate 
their  contributions  to  p ,x  back  to  tl  8  before  the  field-solve: 

Time  level  --3  4  5  6  7  8  9  10  11 


After  the  active  particles  have  been  pushed  to  tl  7  (and  before  their  x- 
contributions  have  been  accumulated),  they  are  moved  if  their  (x,v)  so  dictate 
into  new  blocks.  The  re-distribution  moves  particles  only  into  blocks  that 
will  be  "pre-pushed"  on  the  current  step.  Then  the  pre-push  to  tl  8  (or 
beyond)  is  performed.  Finally,  p  and  x  iare  nterpolated  to  tl  8  (for  all 
blocks,  both  active  and  inactive,  pre-pushed  beyond  8),  and  the  field  equation 
is  solved  for  E8. 

We  also  must  specify  what  At  to  use  in  the  definition  of  the  contribution 
to  x  from  a  block;  assuming  linear  interpolation  of  both  p  and  x.  it  is  correct 
to  use  the  At  associated  with  that  block  (roughly,  also  asssociated  with  a 
region  of  (phase)  space  or  a  grid  spacing).  Thus  the  At  that  goes  into  the 
formula  for  x  is  not  a  constant. 

The  algorithm  may  subsume  explicit  ion- subcycling  schemes;  ions  would 
normally  be  processed  only  at  the  longest  interval.  However,  initialization  is 
simplified  if  all  particles  are  placed  in  blocks  el  or  il  at  the  start,  to 
avoid  referencing  a  negative  time-level.  During  the  first  few  steps  some  delay 
in  "promoting”  particles  into  higher  number  iw  ..'jst  be  imposed,  so  that 
(e.g.)  blocks  i4-7  are  uniformly  populated.  After  a  few  steps,  blocks  il-3 
will  be  empty. 

The  Revised  D1  Particle -Advance  Algorithm 

We  seek  a  variant  of  the  D1  scheme  with  x,v  defined  at  the  same  (integer) 
time  level,  to  facilitate  changing  the  timestep  size.  We  start  with  the  scheme 
as  it  is  usually  written.  The  "final  push”  is: 

an-i  “  |  t^n-a  +  (q/n>)En(Sn)] 

vn-i/2  “  V1/2  +  (At/2)(q/m)En(xn) 

Xyj  “  ^  +  (At2/2)(q/m)En(xn)  . 

Then,  the  "pre-push"  is: 

Vn+l/2  “  vn-l/2  +  (ht/2)an_1 
*n+l  -  *n  +  At^n+l/2  • 
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We  move  the  computation  of  xn+1  to  the  beginning  of  the  "final  push",  where  it 
becomes:  xn  -  xn_1  +  Atvn_1/2.  We  then  relabel  vn+1/2,  calling  it  vn  (it's 
actually  centered  at  tl  n,  formally,  so  this  is  a  notational  improvement) . 

Then,  we  can  write: 

vn  -  Vn-l/2  +  (^t/2)an_1 

“  vn-i  +  (At/2)an.1  +  (At/2)  (q/m)En(xn)  . 

The  Algorithm  in  its  Entirety 

We  enter  a  timestep  with  the  particle  data  x^.j,  vn_lt  and  an_2, 
and  with  on  the  mesh.  Strictly  speaking,  we  should  write  a  trivial 
generalization  of  the  following,  with  incoming  x  defined  at  time  level  n-j , 
etc.,  but  we  write  the  algorithm  as  if  j  were  unity  for  clarity.  In  the 
following,  time- subscripted  quantities  are  stored  in  the  particle  arrays,  while 
unsubscripted  quantities  are  used  only  as  scratch  within  the  particle  loops. 
Explicitly,  the  algorithm  is: 

BEGIN  THE  "FINAL- PUSH”  LOOP  OVER  BLOCKS  AND  PARTICLES: 

0)  At  -  At  (current  block) 

I)  x  -  Vi  +  AtVn-l 
^old  ”  Si -2 

3)  a  -  (q/m)En(x)  (interpolation  of  field  from  mesh) 

4)  S-i  -  \  (a  +  S-21 

5)  vn  ”  Vn-1  +  (At/2)an.1  +  (At/2)a 

6)  =»  x  +  (Atz/2)a 

7)  Enforce  particle  b.c.'s;  reflect  or  absorb  or  shift  by  one  period. 

8)  If  particle  has  moved  to  a  point  in  phase  space  where  At(xn.vn^  —  Atbiock/2' 
set  a^  -  |  [a  +  an_x] ,  and  set  "new  block"  flag. 

9)  If  particle  has  moved  to  a  point  in  phase  space  where  At(xn,vn)  >  2Atblock , 
set  an_1  -  aold,  and  set  "new  block"  flag. 

At  this  point  we  EXIT  THE  "FINAL- PUSH"  LOOP. 

We  then  treat  the  special  cases: 

10)  Sort  flagged  particles  into  new  blocks. 

II)  Inject  any  new  particles  by  adding  them  to  the  appropriate  blocks. 
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At  this  point  we  BEGIN  THE  "PRE-PUSH"  LOOP  OVER  BLOCKS  AND  PARTICLES: 


12)  Update  p,  x  arrays  associated  with  the  current  block  at  tl  n,  using  data 
from  the  individual  particles  which  change  blocks. 

13)  x  -  +  Atblock  vn 

14)  Using  x,  compute  p,  x  arrays  associated  with  the  current  block  at  tl  n+1 . 
At  this  point  we  EXIT  THE  "PRE-PUSH"  LOOP.  We  then: 

15)  Interpolate  p ,x  from  all  necessary  blocks  to  tl  n+1. 

16)  Execute  the  field-solve  to  obtain  En+1. 
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